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ABSTRACT 

This  report  describes  two  computer  programs  that  use  normal  mode 
theory  to  predict  acoustic  propagation  and  discusses  the  use  of  these 
programs  in  analyzing  explosive  data  taken  in  a  series  of  shallow  water 
tests.  Program  S  144 1  calculates  normal  mode  propagation  in  a  medium  in 
which  the  stratification  is  constant  with  range;  program  S  1 548  computes 
normal  mode  propagation  when  the  horizontal  stratification  is  assumed  to 
vary  slowly  with  range.  The  results  of  the  tests  indicate  that  the  introduc¬ 
tion  of  horizontal  stratification  in  program  S  1548  provides  an  improved 
comparison  with  measured  data.  Physical  explanations  are  provided  for 
the  previously  observed  minimum  in  propagation  loss  at  about  125  Hz 
and  for  the  nonlinear  relationship,  at  certain  frequencies,  between  prop¬ 
agation  loss  and  size  of  the  negative  sound  gradient  observed.  In  addition, 
bottom  loss  was  determined  in  the  BIFI  (Block  Island  Fishers  Island)  range 
for  a  wide  range  of  frequencies  and  thermal  conditions. 
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COMPUTER  PROGRAMS  TO  CALCULATE  NORMAL  MODE  PROPAGATION 
AND  APPLICATIONS  TO  ANALYSIS  OF  EXPLOSIVE 
SOUND  DATA  IN  THE  BIFI  RANGE 

INTRODUCTION 

Past  studies  have  shown  that  normal  mode  theory  can  be  useful  in  the  prediction 
of  acoustic  propagation.  This  report  describes  two  FORTRAN  V  computer  programs, 
S1441  and  S1548,  that  can  be  used  for  this  purpose  and  describes  the  use  of  these 
programs  in  the  analysis  of  explosive  data  taken  in  a  series  of  shallow  water  acous¬ 
tic  tests. 

Program  S1441  (discussed  in  section  1  of  this  report)  is  an  extension  of  a  pro¬ 
gram,  written  for  the  Navy  Underwater  Sound  Laboratory*  by  A.  D.  Little,  Inc. 
(references  1  and  2),  that  considers  normal  mode  propagation  in  a  medium  in  which 
the  stratification  is  constant  with  range.  In  the  A,  D.  Little  program,  the  ampli¬ 
tude  distribution  of  an  acoustic  signal  as  a  function  of  depth  is  determined  for  a 
given  mode  by  means  of  the  numerical  solution  of  the  acoustic  wave  equation  for 
given  boundary  conditions.  Program  S1441  extends  the  A.  D.  Little  program  by 
calculating  and  producing  calcomp  plots  of  the  following  values  for  any  mode, 
frequency,  and  velocity  profile: 

a.  Amplitude  as  a  function  of  depth,  and  the  ray  equivalent  of  any  mode 

b.  Group  velocity,  phase  velocity,  excitation  pressure,  and  the  angle  of 
incidence  of  sound  waves  striking  the  boundaries 

c.  Propagation  loss  as  a  function  of  range. 

Program  SI 548  (discussed  in  section  2  of  this  report)  uses  normal  mode  theory 
to  predict  acoustic  propagation,  in  a  medium  whose  velocity  profile  varies  slowly 
with  distance  from  the  acoustic  source,  over  an  ocean  bottom  whose  depth  and 
acoustic  impedance  change  slowly  with  range.  It  calculates  and  produces  calcomp 
plots  of  the  following  values  for  given  modes  at  any  frequency: 


'Now  the  New  London  Laboratory  of  the  Naval  Underwater  Systems  Center  (NUSC). 
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a.  Amplitude  as  a  function  of  depth,  and  the  ray  equivalent  of  any  mode  at 
given  distances  from  a  source 

b.  Propagation  loss  as  a  function  of  range. 

The  two  programs  were  used  in  the  analysis  of  data  taken  in  a  series  of  acous¬ 
tic  tests  conducted  in  Block  Island  Sound  between  August  1967  and  October  1968. 
In  these  tests,  referred  to  as  "Experiment  2,f  in  reference  3,  propagation  loss  was 
measured  under  a  wide  range  of  thermal  conditions,  using  explosives  as  sound 
sources.  In  section  3  of  this  report,  these  propagation  loss  measurements  are  com¬ 
pared  with  the  values  predicted  by  normal  mode  theory. 

This  report  is  based  on  material  contained  in 

•  W,  G.  Kanabis,  "A  Computer  Program  to  Calculate  Normal  Mode  Propa¬ 
gation  Over  a  Flat  Homogeneous  Ocean  Bottom,11  NUSL  Technical  Memorandum 
No.  2211-296-69,  13  October  1969. 

•  W.  G.  Kanabis,  "A  Computer  Program  to  Calculate  Normal  Mode  Propa¬ 
gation  in  a  Medium  in  Which  Stratification  Is  a  Function  of  Position, 11  NUSL 
Technical  Memorandum  No.  2211-11-71,  14  January  1971. 

•  W.  G.  Kanabis,  "Propagation  of  Explosive  Sound  in  the  BIFI  Range," 
NUSL  Technical  Memorandum  No.  2211-311-70,  9  November  1970. 
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Section  1 

A  COMPUTER  PROGRAM  (S 1 44 1 )  TO  CALCULATE  NORMAL  MODE 
PROPAGATION  OVER  A  FLAT,  HOMOGENEOUS  BOTTOM 

NORMAL  MODE  THEORY 


Normal  mode  theory  is  based  on  the  assumption  that  at  large  distances  from  a 
source  the  main  part  of  the  sound  field  consists  of  standing  waves  formed  by  energy 
striking  the  boundaries  in  certain  preferred  directions. 


There  are  three  methods  by  which  the  preferred  directions  and  the  amplitude 
distribution  of  the  standing  waves  may  be  determined.  First,  one  may  find  the 
solution  of  the  wave  equation  in  closed  form: 


V24>  = 


1  d2<t> 


2  a  2 
c  at 


0) 


where  $  is  a  displacement  or  a  velocity  potential  and  c  is  the  veloci ty  of  sound 
in  the  medium  considered.  This  is  done  in  reference  4  by  integrating  equation  (1) 
subject  to  the  boundary  conditions  in  the  complex  plane  and  approximating  the 
normal  mode  solution  for  the  standing  waves  by  the  evaluation  obtained  from  the 
residues  in  the  integration. 


Second,  one  may  obtain  the  solution  of  equation  (1)  by  direct  numerical  inte¬ 
gration,  as  done  by  A.  D.  Little,  Inc. ,  by  means  of  the  computer  program  described 
in  references  1  and  2. 


Third,  as  shown  in  reference  5,  one  may  consider  as  the  descriptions  of  the 
modes  those  ray  paths  that  undergo  total  reflection  at  the  boundaries  and  whose 
successive  upgoing  and  downgoing  rays  are  in  phase. 


SOLUTION  OF  THE  WAVE  EQUATION  BY  DIRECT 
NUMERICAL  INTEGRATION 

The  basis  of  this  program  is  the  rapid  evaluation  afforded  by  the  use  of  a  high¬ 
speed  computer. 

When  equation  (1)  is  solved  for  4>,  the  incremental  pressure  p  is  found  by 
definition 
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d 2<t> 

P  P  2  ' 

dt 


(2) 


where  4>  is  the  displacement  potential  and  p  is  the  material  density  at  the  point 
considered . 

We  find  that  equation  (1)  is  separable  in  terms  of  range  r  and  depth  z,  so 
that  the  pressure  amplitude  can  be  written  as 


Pa  =  F(r)  •  u(z) 

and  a  differential  equation  in  terms  of  u(z)  may  be  obtained: 


(3) 


(4) 


where 

u  is  the  radial  frequency, 
c  is  the  sound  velocity  (a  function  of  z),* 

u  =  u(z)  is  the  pressure  amplitude  distribution  as  a  function  of  depth,  and 

kr  is  the  horizontal  wave  number 

—  sin  0  ,  (5) 

c 

in  which  0  is  measured  relative  to  the  normal  to  the  ocean  bottom  and  c  and  0 
are  measured  at  the  same  depth. 

The  physical  picture  presented  by  equation  (3)  is  of  a  standing  wave  that  has 
a  particular  pressure  amplitude  distribution  as  a  function  of  depth,  u(z),  and  that 
travels  unchanged  in  shape  as  it  progresses  in  the  r  direction. 


Equation  (4)  may  be  written  as 

d2u 

dz2 


+  f(z)u  =  0  , 


(6) 


*The  z  dependence  will  be  dropped  henceforth  in  the  notation  c. 
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where 


(7) 


In  solving  equation  (4),  we  must  impose  two  constraints  on  any  solution:  the 
boundary  conditions  at  the  interface  between  the  water  and  the  bottom  material 
and  at  the  water  surface.  First,  because  there  must  be  continuity  of  pressure  and 
particle  velocity  in  the  vertical  direction  at  the  bottom  interface,  it  is  necessary 
that 


_L 

u  i  dz 


(8a) 


where  the  terms  with  subscript  b  refer  to  quantities  in  the  bottom  material  while 
those  with  subscript  1  refer  to  the  water  side  of  the  interface.  Second,  since  we 
assume  a  pressure -re  I  ease  surface  at  the  air-water  boundary,  then 

u  =  0  (8b) 


at  this  surface . 

Also,  we  assume  that  only  "unattenuated"  modes  compose  the  sound  field. 
These  modes,  by  definition,  involve  propagation  with  energy  that  strikes  the  bot¬ 
tom  at  angles  larger  than  the  critical  angle  so  that  "total  reflection"  occurs. 
Since  the  critical  angle  0C  measured  relative  to  the  normal  to  the  interface  is 
given  by 

sin  9C  =  — -  , 
b 

then 

sin  9  >  — * — 

Cb 

must  hold  for  all  energy  striking  the  bottom.  (Naturally,  must  be  greater  than 
to  ensure  the  existence  of  a  critical  angle.) 
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Since 


k 


r 


t 


using  equation  (7)  and  the  above  inequality,  we  obtain 


2  2 
u  to 


I'Wl  <  ^2  -  ^ 


C1  Cb 


(8c) 


at  the  bottom  interface. 

Upon  examining  equation  (6),  we  can  see  that  a  solution  of  the  equation  is 


u(z)  =  B(z)e 


iVffe) 


for  given  values  of  c,  to,  and  kr.  B  is  a  constant.  If  f(z)  is  positive,  u(z)  is  in 
the  form  of  a  sinusoid .  This  form  of  solution  is  obtained  for  the  interference  pattern 
between  upgoing  and  downgoing  waves  in  the  water.  If  f(z)  is  negative,  u(z)  is 
in  the  form  of  an  exponential .  This  is  a  consequence  of  Sturm^  comparison  theorem 
(reference  6). 


The  value  f(z)  is  negative  everywhere  in  the  bottom  and  in  the  water  at 
depths  that  correspond  to  shadow  zones  caused  by  the  vertfcxing  of  rays  that  form  a 
mode.  Both  distributions  are  a  result  of  the  condition  of  continuity  of  pressure  in 
the  medium.  Thus,  when  there  is  total  reflection  at  a  level  (either  by  vertexing 
or  reflection  from  a  boundary)  and  the  pressure  is  finite  at  that  level,  then  at 
adjacent  levels  there  is  a  decay  (whose  rate  is  determined  by  the  boundary  con¬ 
ditions)  but  not  a  discontinuous  step  to  zero  pressure.  Examples  of  a  distribution 
involving  both  exponentials  are  shown  in  figure  1/  which  presents  the  amplitude 
distribution  as  a  function  of  depth,  and  in  figure  2,  which  presents  the  ray  equiva¬ 
lent  of  the  mode.  This  concept  of  the  ray  equivalent  is  discussed  in  more  detail  in 
the  next  subsection  of  this  report.  In  figure  2,  it  is  seen  that  the  ray  equivalent  of 
the  normal  mode  vertexes  at  a  depth  of  48  ft;  therefore,  in  figure  1,  the  pressure 
amplitude  distribution  is  in  the  form  of  an  exponential  between  the  surface  of  the 
ocean  and  a  depth  of  48  ft.  If  can  also  be  seen  i  n  figure  1  that  the  pressure  amplitude 
distribution  in  the  bottom  is  in  the  form  of  an  exponential. 

Figures  1  through  16  oppeor  ot  the  end  of  this  section,  pages  23  through  33. 
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To  integrate  equation  (6)  numerically,  formulas  relating  un  +  ^  and  its  deriv¬ 
ative  Up  +  ]  with  the  quantities  un  and  u^  must  be  obtained  (the  subscripts  sig¬ 
nify  the  depth  at  which  the  quantities  are  calculated).  This  is  done  by  writing  a 
Taylor  series  for  un+]  and  Un+1  and  retaining  al I  terms  of  the  third  order  or 
less.  The  details  of  this  procedure  are  given  in  references  1  and  2;  the  results  are 
summarized  in  equations  (9)  and  (10): 


u 


n  +  1 


un  +  hun 


(9) 


u  ,  1 
n  +  1 


+  f. 


n  +  1 


n  +  1 


(10) 


where 

h  is  the  increment  between  level  n  and  level  n  +  1  and 

fn  and  fp  +  ^  are  the  values  of  f(z)  at  n  and  n+1,  respectively. 

Equations  (9)  and  (10)  are  recursion  relationships  that  permit  calculation  of  u 
and  u1  at  all  levels  in  the  water  if  the  given  values  for  u^  and  u*i  are  the  values 
of  u  and  u1  just  above  the  bottom. 

Since  we  are  interested  in  normalized  values  of  u  over  the  water  column,  we 
can  select  u^  to  be  any  arbitrary  value  (u ^  =  1  is  convenient) .  For  an  arbi trary 
value  of  the  horizontal  wave  number  kr  that  is  restricted  by  equation  (8c),  we  can 
evaluate  u\  by  equation  (8a) .  Then  we  can  determine  u  for  all  levels  by  using 
equations  (9)  and  (10)  repeatedly.  If  the  value  of  kf  corresponds  to  a  mode,  equa¬ 
tion  (8b)  will  be  satisfied  at  the  surface  of  the  water.  Each  mode  has,  at  most,  one 
such  solution  for  a  given  frequency.  There  is  a  low  cut-off  frequency  for  each  mode, 
sothatat  frequencies  below  the  cut-off  frequency,  equation  (8b)  can  not  be  satisfied. 

After  finding  the  amplitude  distribution  of  a  mode,  we  must  define  the  mode 
number.  For  a  finite  frequency,  the  mode  number  is  equal  to  the  number  of  nodes 
in  the  amplitude  distribution.  Thus,  the  first  mode  has  a  node  only  at  the  surface. 

A  representation  of  the  amplitude  distribution  of  the  first  mode  is  shown  in  figure  1 . 
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THE  RAY  EQUIVALENT 


Corresponding  to  the  definition  of  a  mode  discussed  above  is  a  more  physical 
approach  in  which  the  ray  equivalent  of  the  solution  to  the  wave  equation  is  con¬ 
sidered.  Additional  information  concerning  this  approach  is  contained  in  references 
4,  5,  7,  8,  and  9. 

For  simplicity,  let  us  consider  a  two-layer  medium  of  constant  water  depth  H, 
density  p-j,  and  sound  velocity  c-j  lying  over  an  infinite  bottom  of  density  P2 
and  sound  velocity  C2,  as  shown  in  figure  3.  At  large  ranges  from  a  point  source, 
we  may  consider  sound  to  be  propagated  by  plane  waves.  It  is  clear  from  figure  3 

that  for  certain  waves  whose  direction  is  defined  by  an  angle  0,  there  will  be 

constructive  interference  between  that  wave  and  a  plane  wave  undergoing  one 
more  bottom  and  surface  reflection.  For  constructive  interference  to  exist,  the 
phase  difference  between  points  A  and  B  of  figure  3  must  be  2(n  -  1)tt  degrees;  i.e., 
the  phase  difference  must  satisfy  the  equation 

(~  H  Q  +  H  COS  20]  -  e  -  ir  =  2(n  -  l)ir 

An  [cos  0  cos  0  J 


or 


[2Hcos  0]  -  e  -tt  =  2(n-l)ir  ,  (11) 

where  Xn  is  the  wavelength  of  the  preferred  mode,  e  is  the  phase  change  under¬ 
gone  by  a  plane  wave  upon  bottom  reflection,  n  is  the  mode  number,  and  the  phase 
change  upon  reflection  from  the  water  surface  is  assumed  to  be  -tt  .  If  the  sound 
velocity  in  the  water  layer  varies  with  depth,  the  first  term  of  equation  (1 1)  would 
be  different  from  that  given  above;  however,  the  discussion  below  applies  to  either 
case . 

If,  for  a  given  wavelength  Xn  and  angle  0,  there  is  constructive  interference 
between  plane  waves  suffering  different  numbers  of  bottom  and  surface  reflections, 
then  propagation  consists  of  a  series  of  upgoing  and  downgoing  waves,  as  shown  in 
figure  4.  The  left-hand  term  of  equation  (11)  is  the  phase  change,  2A,  undergone 
in  the  z  direction  when  a  ray  makes  a  surface-bottom-surface  cycle.  For  finite 
frequencies,  tt  >  €  >  0  as  tt/2>0>0c.  Therefore,  the  phase  change.  A,  under¬ 
gone  in  the  z  direction  over  the  water  depth  is  limited  by 

A  <  mr  (12) 
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The  pressure  is  zero  at  the  surface,  the  phase  change  upon  reflection  is  -a,  and 
the  direction  of  propagation  is  reversed  upon  reflection  from  the  surface.  There¬ 
fore,  the  sound  field  in  the  vertical  direction  is  the  sum  of  two  sine  waves  repre¬ 
senting  the  upgoing  and  downgoing  waves  in  the  z  direction.  These  waves  are 
shown  for  the  first  two  modes  in  figure  5.  Because  of  equation  (12),  the  number  of 
nodes  in  this  amplitude  distribution  is  equal  to  the  mode  number.  Thus,  there  is  a 
correspondence  between  the  definitions  of  mode  number  in  the  solution  of  the  wave 
equation  and  the  ray  equivalent  solution. 

Now  let  us  consider  the  procedure  actually  used  in  the  computer  programs  to 
determine  the  ray  equivalent.  When  the  wave  equation  is  solved  numerically, 
values  of  f(z)  are  obtained;  f(z)  is  given  in  equation  (7)  by 

f(z)  —  ,  (13) 

c 

where  kr  is  given  in  equation  (5)  by 

kr  =  sin  0  .  (14) 

Therefore,  given  positive  f(z),  one  can  determine,  from  equations  (13)  and 
(14),  the  cosine  of  the  angle  of  inclination  of  the  equivalent  ray  as  a  function  of  z : 

cos  9  =  77  ^  f(z)  •  (15) 

U 


If  the  ray  between  two  points  z^  and  Z2  is  continuous,  then  Az  =  z^  -  Z2 
may  be  given  in  terms  of  the  horizontal  distance  AR  =  -  R2  and  one  particular 

value  of  the  tan  8  over  the  path.  This  relationship  is 

AR  n  c  sin  0  c 

-7 —  =  tan  0  =  - ;r  - - . 

Az  c  cos  0  c  cos  0 


where  c  ,  the  vertex  velocity,  equals  c/sin  0. 


If  the  value  of  8  does  not  vary  appreciably  between  z^  and  Z2,  then  we 
can  approximate  tan  0  by 

tan  0i  +  tan  89 
tan  0= - - -  . 
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where  0-|  and  02  are  taken  at  z  =  z-|  and  respectively.  Thus  (reference  10), 


AR  _  1  c1  +c2 
Az  cv  cos  0-j  +  cos  ©2 


(16) 


where 

cy  is  the  vertex  velocity, 

c-|,  C2  are  the  sound  velocities  at  z-|  and  Z2,  respectively,  and 

©  1 '  @2  are  the  angles  relative  to  the  normal  ofthe  ray  at  z^  and  Z2,  respectively. 

Thus,  given  values  of  f(z),  one  can  construct  a  ray  equivalent.  If  vertexing 
takes  place,  the  depth  zv  at  which  it  occurs  is  the  level  whose  value  of  sound 
velocity  equals  cy. 

The  ray  equivalent  has  a  counterpart  in  the  analysis  of  sound  propagation  using 
ray  theory.  In  ray  theory,  rays  corresponding  to  a  continuum  of  angles  of  propaga¬ 
tion  are  summed  at  the  receiver.  Most  groups  of  rays  effectively  cancel  each  other 
and  contribute  little  to  the  sound  field.  This  leaves  a  discrete  number  of  rays  that 
sum  constructively  to  form  the  sound  field.  These  rays  correspond  to  the  ray  equiv¬ 
alent  of  modes  that  compose  the  sound  field.  According  to  ray  theory,  the  rays  that 
contribute  to  the  sound  field  are  determined  by  interference  effects,  and  hence,  as 
in  normal  mode  theory,  the  geometry  of  the  dominant  rays  changes  as  a  function  of 
frequency.  Also  as  in  normal  mode  theory,  the  sum  of  these  rays  is  influenced  by 
the  source-receiver  geometry.  However,  there  are  two  factors  that  limit  the  corre¬ 
spondence  between  normal  mode  and  ray  theory.  First,  in  ray  calculations  the  sound 
field  is  often  obtained  by  simply  adding  all  rays  without  regard  to  their  phase. 
Second,  the  geometrical  approximation  limits  the  validity  of  ray  tracing.  This 
limitation  can  produce  differences  not  only  in  the  theoretical  pressure  field  but 
also  in  the  ray  representations  in  ray  tracing  and  normal  mode  methods. 

If  the  velocity  profile  assumed  in  the  calculations  contains  only  a  monotonic 
variation  in  sound  velocity  with  depth  (as,  for  example,  in  figure  6)  or  describes  a 
simple  sound  channel,  a  relatively  simple  ray  equivalent  obtains  (as  in  figure  2), 
However,  let  us  consider  the  ray  equivalent  associated  with  the  more  complicated 
velocity  profile  shown  in  figure  7.  This  profile  describes  two  sound  channels  at 
depths  of  about  60  and  120  ft,  respectively.  The  associated  ray  equivalent  of 
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mode  1  at  a  frequency  of  200  Hz  is  shown  in  figure  8.  It  can  be  seen  that  the  ray 
equivalent  consists  of  two  traveling  waves,  one  for  each  of  the  two  sound  channels 
in  the  water  column.  This  phenomenon  of  wave  "leakage"  through  a  layer  is  des¬ 
cribed  in  reference  11. 

PHASE  VELOCITY  AND  GROUP  VELOCITY 

The  phase  velocity  Vp  is  given  by  the  relationship 

V  =  c 
P  sin0  ' 


07) 


where  0  is  the  direction  of  propagation  of  a  plane  wave  where  the  sound  velocity 
has  a  value  c.  Equation  (17)  can  also  be  written 


where  cy  is  the  vertex  velocity. 


08) 


The  group  velocity  Vg  can  be  considered  from  two  points  of  view.  First,  V 
may  be  considered  as  a  measure  of  the  speed  of  propagation  in  the  horizontal  direc 
tionof  a  number  of  frequencies  in  a  band  Au  centered  about  u.  Vg  may  be  given 
by  (reference  12) 


(19) 


or  by  (reference  11) 


V 


9 


Vr 


c  c 


+ 


(20) 


where  H  is  the  water  depth. 

It  can  be  seen  that  this  approach  to  the  calculation  of  group  velocities  in¬ 
volves  the  calculation  of  derivatives,  which  is  not  desirable  in  a  computer  program 
since  it  produces  inaccuracies  and  makes  it  necessary  to  obtain  an  unnecessarily 
large  number  of  values  of  phase  velocity  as  a  function  of  frequency. 


11 


I  CO 


TR  4319 


Second,  Tolstoy  (reference  13)  has  used  a  general  theorem  by  Biot  (reference  14) 
to  show  the  equivalence  of  Vg  in  equations  (19)  and  (20)  to  the  rate  of  energy 
transport  in  the  horizontal  direction.  The  group  velocity  is  given  by 


V 


g 


(21) 


where 


v 


<j 


f. 

+ 

I. 


dz 


JL  4>2  dz 

C 


(22) 


(23) 


in  which  p  and  c  are,  respectively,  the  density  and  sound  velocity  at  z,  and 
<f>  is  given  in  the  equation 

<|>  =  <{>  (z)  e  (24) 

in  which  <t>  is  the  displacement  potential  in  equation  (1). 


Since,  by  equation  (2), 


P  =  ~P 


d2$ 

at2 


and  u  is  the  value  of  pressure  (normalized  to  maximum  amplitude)  as  a  function 
of  z;  then,  by  equations  (2)  and  (24), 

uapip  ,  (25) 

where  <P  is  the  normalized  value  of  cp.  Thus,  by  equation  (25)  we  can  obtain  <p 
once  u  is  known  since  we  are  interested  only  in  the  normalized  values  of  u  and 
<p  for  given  u. 

Given  p-]  in  the  water,  in  the  bottom,  and  u  normalized  to  maximum 
pressure  amplitude;  to  obtain  ip,  the  normalization  of  <p,  we  first  obtain 
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„  .  “W. 

1/uu(z)  Pi 
_  u(z) 

b,  %(z)  pb  ' 


(25a) 


where  1  and  b  signify  water  and  bottom,  respectively,  and  u  signifies  unnor¬ 
malized.  Since  the  maximum  value  of  u  lies  in  the  water  and  p^  is  greater  than 
p^  and  because  u(z)  is  normalized  with  respect  tomaximum  amplitude,  wemultiply 
both  expressions  in  equation  (25a)  by  p^  to  obtain  ip  normalized  with  respect  to 
maximum  amplitude,  so  that 


V»i  =  u(z) 

%  =  u(z)  P]/ Pb 


(25b) 


The  normalized  value  for  <f>  is  represented  by  ip  and  may  be  used  in  place  of  <f> 
in  equations  (21)  through  (23). 

EXCITATION  PRESSURE  AND  PROPAGATION  LOSS 


The  sound  field  produced  by  a  simple  harmonic  source  in  a  two-layered  half¬ 
space  (figure  9)  with  a  free  surface  at  z  =  0  and  the  boundary  between  two  fluids 
at  z  =  H  is  given  by  the  solution  of  equation  (1).  The  solution  is  given  in  refer¬ 
ence  12  by 


1  1 
Vr~  u2 


m 


p  -iUmr-ut-n/4) 

rme 


(26) 


where 

r  is  the  horizontal  range, 
cj  is  the  radial  frequency, 
m  is  the  mode  number, 

is  the  horizontal  wave  number  kf  for  mode  m,  and 

^m  ~~  Pm  ps  ^m^zo^  '  (^7) 
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in  which 

ps  is  the  water  density  at  the  source, 

<^m(z)  is  the  normalized  displacement  potential,  a  function  of  depth, 
zQ  is  the  source  depth, 
z  is  the  receiver  depth,  and 
m  is  the  mode  number. 

The  notation  pm  is  the  excitation  function  given  by 

Pm  =  2ir  (PocoS)1/2  - T==~  '  (28) 

■'m  V*m 

where 

pQ  is  the  water  density  at  the  source, 

cQ  is  sound  velocity  at  the  source, 

S  is  the  power  output  of  an  omnidirectional  source,* 

k  is  the  horizontal  wave  number,  and 
m 


4-00 


vm 


L 


2 

PVV  dz 


m 


(29) 


where  is  the  normalized  displacement  potential. 


If  the  source  produces  a  unit  sound  pressure  level ,  then  the  following  relation¬ 
ship  (reference  11)  must  be  satisfied: 


4tt  I 


Sco  Pc 
v  2ir 


1/2 


=  1 


(30) 


Substituting  equation  (30)  into  equation  (28),  we  obtain  for  the  excitation  pressure 


*This  quantity  is  represented  by  II  in  reference  12. 
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P 


m 


,0  \  1/2 

P0(2*) 

2vmV\7 


(31) 


The  excitation  pressure  is  the  sound  pressure  amplitude  produced  by  a  source  that 
generates  a  unit  pressure  level  at  unit  range  when  both  source  and  receiver  are  at 
antinodes.  It  is  essentially  a  measure  of  the  source  level  of  mode  m  for  a  unit  source. 


From  equation  (26)  we  determine  the  pressure  amplitude  characteristics  of  the 
d  field,  and  this  amplitude,  p  ,  is  given  by 

1/2 

P»=^{[SPm'“{'mr-T)]  +  [SPmsin  (Vi)]  )  '  (32) 


Since  pQ  is  the  sound  pressure  amplitude  at  range  r  from  a  generator  with 
unit  source  level,  the  value  of  propagation  loss  Lf  at  range  r  for  given  depths 
of  source  and  receiver  is 


Lr  =  -20  log  pQ 


(33) 


It  can  be  seen  from  equations  (27)  and  (32)  that,  once  is  known,  one  can 
easily  determine  the  effect  of  the  source  and  receiver  depths  on  the  sound  field  at 
a  given  range.  If  the  source  depth,  zQ,  is  such  that  ^m(z0)  is  a  node,  the  mode 
m  will  be  suppresed  in  the  sound  field;  conversely,  if  zQ  is  such  that  v?m(z0)  is 
an  antinode,  the  sound  field  of  mode  m  will  be  greater  than  it  would  be  at  depths 
for  which  <^m(z)  is  less  than  <pm(zQ).  The  same  relationships  apply  to  the  effect 
of  the  receiver  depth  upon  the  sound  field. 

It  can  also  be  seen  from  equations  (27),  (28),  and  (29)  that  the  pressure  ampli¬ 
tude  does  not  depend  upon  the  normalization  of  cp(z) . 

For  small  attenuations  of  individual  modes,  equation  (32)  may  be  rewritten  to 
include,  for  mode  m,  losses  at  the  boundaries  as  a  function  of  range  and  losses 
caused  by  absorption  of  sound  energy  in  the  water,  so  that 


Pa 


EP  10(~Dmr/2°-ar/20) 

m  m 


+ 


Ep  ,0<-Dmr/»-°'/20) 

m  m 


sin  Ik  r 
m 


2 

1/2 


(34) 
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where  Dm*  is  a  measure  of  bottom  loss  in  decibels  of  loss  per  unit  increment  of 
range  and  a  is  the  attenuation  coefficient  given  in  reference  15  by 


a 


dB 

kyd 


0.1  f2  40  f2 
1  +  f2  4100  +  f2 


where  f  is  the  frequency  in  kilohertz. 

DESCRIPTION  OF  PROGRAM  S 1441 

Program  S1441  uses  normal  mode  theory  to  predict  acoustic  propagation  over  a 
flat,  homogeneous  ocean  bottom.  The  program  may  be  used  with  any  of  three  options, 
each  of  which  provides  different  information  about  the  sound  field  in  a  medium  for 
a  given  frequency,  velocity  profile,  and  mode  number.  The  three  options  provide, 
as  output,  the  following  calcomp  plots. 

•Option  A  produces  two  types  of  plots  for  each  mode  analyzed.  One  type 
gives  pressure  normalized  to  the  maximum  amplitude  as  a  function  of  depth, 
and  the  second  gives  the  ray  equivalent  of  the  mode. 

•Option  B  produces  two  types  of  plots  for  each  mode.  One  type  gives  three 
quantities:  phase  velocity,  group  velocity,  and  excitation  pressure  as  a 
function  of  frequency,  and  the  second  gives  the  angles  of  incidence  of 
energy  at  the  two  boundaries  as  a  function  of  frequency. 

•Option  C  produces  a  plot  of  propagation  loss  versus  range  for  any  combi  nation 
of  modes.  Plots  can  be  produced  for  any  source  depth,  receiver  depth,  and 
frequency. 

Each  option  provides  two  plots  of  the  velocity  profile:  the  first  plot  (figure  6) 
shows  the  sound  velocity  in  both  the  water  and  the  bottom;  the  second  (figure  10) 
shows  the  sound  velocity  in  the  water  in  greater  detail  through  the  use  of  an  ex¬ 
panded  velocity  scale. 

The  format  of  the  input  data  is  presented  in  table  1  .  Card  group  2  in  the  table 
shows  the  routine  for  selecting  option  A,  B,  or  C.  The  mechanics  of  the  individual 
options  are  described  below. 


* 


must  be  specified  by  the  user  of  the  program. 
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Table  1.  Input  Data  for  Program  SI 441 


Card 

Group 

Amount 
of  Cards 

Format 

Input 

Parameter 

Columns 

Data 

1 

1 

20A4 

1-80 

Heading  on  each  page  of  output. 

2 

1 

110 

NUMV 

1-10 

Number  of  velocities  in  profile. 

F10.0 

VEL  1 

11-20 

Velocity  at  origin  of  calcomp 

plot  of  velocity  profile  (ft/sec) 

4110 

IVOP 

21-30 

If  IRP  =  0, 

when  IVOP  =  0,  option  A 
when  IVOP  =  1,  option  B 

IRP 

31-40 

When  IVOP  =  1, 

i f  IRP  =  1 ,  option  C 

IVRP 

41-50 

When  IVRP  =  1,  velocity  profile 
not  plotted 

When  IVRP  =0,  profile  plotted. 

1 

r 

IEX 

51-60 

Changes  increment  of  kr  (hori¬ 

zontal  wave  number)  by  a 
factor  of  10-IEX#  Values  from 

0  to  10. 

3 

1 

5F 

0.3 

ZM 

1-10 

Water  depth  (ft) 

CB 

11-20 

Velocity  of  sound  in  the  bottom 

(ft/  sec) 

RO 

21-30 

Density  of  water  (grams/cm^) 

RB 

31-40 

Density  of  bottom  (grams/ cm3) 

' 

FSC 

41-50 

Maximum  depth  plotted  (200  x 

FSC  in  all  velocity  profiles 
and  in  plots  in  option  A) 

4 

NUMV 

2F 

0.3 

Z(l) 

1-10 

Height  above  bottom  at  which 

cards 

r 

C(l) 

11-20 

sound  velocity  is  C(l). 

Velocity  of  sound  (ft/sec)  at 

Z(l). 

1  =1,  2, 

. . .,  NUMV  in  order  of  increasing  Z 

5 

1 

110 

N 

1-10 

Number  of  intervals  into  which 

depth  is  to  be  subdivided  in 
integration  of  differential 
equations. 
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Table  1  (Cont'd).  Input  Data  for  Program  S1441 


Card 

Group 

Amount 
of  Cards 

Format 

Input 

Parameter 

Columns 

Data 

5 

2F10.3 

UM 

1  1-20 

Maximum  value  of  u(z) 

FQ 

21-30 

Frequency  (Hz) 

6 

i 

4110 

1SFQ 

1-10 

Highest  frequency  to  be  analyzed 

IEFQ 

11-20 

Lowest  frequency  to  be  analyzed 

NMOD 

21-30 

Number  of  modes  analyzed. 

Mode  number  analyzed  =  1, 

2,  ...,  NMOD. 

1 

r 

INCF 

31-40 

Decrement  in  frequency  from 

ISFQ  to  IEFQ.  (In  options  A 
and  C,  ISFQ  =  IEFQ  =  FQ. 
Option  B  range  of  frequencies 
is  selected  by  ISFQ,  IEFQ, 
INCF.) 

7 

i 

4110 

1RST 

1-10 

Shortest  range  (ft)  at  which 

(for 

propagation  loss  versus  range 

option 

will  be  plotted. 

C) 

IREN 

11-20 

Longest  range  (ft)  at  which 

propagation  loss  will  be  plotted. 

1R1C 

21-30 

Increment  in  range  to  be  plotted 

(ft). 

} 

r 

NPS 

31-40 

Number  of  propagation  loss  versus 

range  plots. 

3F10.3 

ZS 

41-50 

Source  depth  (ft) 

1 

ZRC 

51-60 

Receiver  depth  (ft) 

T 

FM1 

61-70 

Increment  of  range  on  the 
calcomp  plot  per  division  in 
nautical  miles 

8 

the 

4F10.5 

DD(J) 

1-40 

Loss  at  boundaries  (dB/ft)  as  a 

(for 

smal  lest 

function  of  range  forthemode  J. 

option 

integer 

J  =  1 ,  2, . . . ,  NMOD  where 

C) 

>  NMOD 

NMOD  is  the  number  of  modes 

4 

analyzed. 
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Table  1  (Cont'd).  Input  Data  for  Program  S1441 


Card 

Group 

Amount 
of  Cards 

Format 

Input 

Parameter 

Columns 

Data 

9 

Card  a 

no 

NMS 

1-10 

Number  of  modes  to  be  summed 

(for 

=  1 

in  the  plot  of  propagation  loss 

option 

versus  range  (NM  <  NMOD) 

C) 

Card(s)  b 

6110 

MDS(J) 

1-60 

J  =  1,  2,  ....  NMS 

=  the 

The  values  of  MDS(J)  are  the 

smallest 

mode  numbers  of  the  modes  to 

integer 

be  summed  in  the  propagation 

>  NMS 

loss  versus  range  plots. 

6 

MDS(J)  <  NMOD 

(Cards 

a  and  b  form  a  set.  There  is  a  set 

of  cards  for  each  plot  in  option  C.) 

OPTION  A 

In  option  A,  first  the  numerical  solution  to  equation  (6)  is  found,  subject  to  the 
boundary  conditions  given  by  equations  (8a),  (8b),  and  (8c);  then  the  ray  equiva¬ 
lent  of  the  solution  is  obtained. 

For  a  given  inputted  velocity  profile  (card  group  4,  table  1),  the  sound  ve loci ty 
is  calculated  at  N  levels  (card  group  5)  equispaced  between  the  surface  and  bot¬ 
tom  by  interpolation. 

The  value  of  the  horizontal  wave  number  kr  is  varied  subject  to  the  restric¬ 
tions  given  in  equation  (8c).  For  each  value  of  kr,  the  value  of  u(z)  is  calculated 
over  the  water  column  of  N  equispaced  levels  by  means  of  equations  (9)  and  (10). 
The  Values  of  u(z)  are  restricted  so  that  u(z)  never  exceeds  UM  (card  group  5). 
In  addition,  mode  m  must  have  m  zero  crossings.  After  u(z)  is  calculated  for  a 
given  kf  value,  kf  is  incremented  by  Akr  so  as  to  obtain  the  smallest  possible 
value  of  u  at  the  surface.  If  the  conditions  necessary  for  the  existence  of  a  given 
mode  can  not  be  met,  the  statement  "No  Mode  Found"  is  printed  out. 

Once  u(z)  has  been  found  for  a  given  mode,  the  ray  equivalent  can  be  found 
by  using  equation  (16).  Sample  calcomp  plots  of  the  amplitude  distribution  normal¬ 
ized  to  maximum  amplitude  and  of  the  ray  equivalent  are  shown  in  figures  1  and  2. 
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Where  there  is  a  large  negative  velocity  gradient  high-frequency  sound  is 
trapped  near  the  ocean  bottom.  Under  these  circumstances,  it  is  difficult  to  obtain 
a  good  approximation  of  the  boundary  condition  at  the  surface,  and  one  has  to  de¬ 
crease  the  value  by  which  kr  is  incremented.  This  is  done  by  using  a  large  value 
(up  to  about  10)  of  IEX  (card  group  2,  table  1).  A  large  value  of  IEX  will  increase 
the  program  time  since  it  decreases  the  increment  of  kr  by  a  factor  of  10“IEX, 
Under  most  circumstances,  however,  an  IEX  value  of  zero  is  adequate. 

OPTION  B 

In  option  B,  group  velocity,  phase  velocity,  excitation  pressure,  and  angle  of 
incidence  at  the  boundaries  relative  to  the  normal  to  the  boundary  are  found  over  any 
frequency  range  (card  group  6,  table  1).  The  group  velocity  for  a  given  frequency  is 
determined  by  equations  (21),  (22),  and  (23);  phase  velocity  by  equation  (18);  ex¬ 
citation  pressure  by  equation  (31);  and  angle  of  incidence  at  the  boundaries  by 
equation  (15) .  If  the  ray  vertexes  before  striking  a  boundary,  the  angle  of  incidence 
is  given  as  90°.  Sample  plots  obtained  from  option  B  are  shown  in  figures  11  and  12. 

OPTION  C 

In  option  C,  propagation  loss  for  a  source  level  at  a  1  —yd  reference  is  obtained 
as  a  function  of  range  for  any  given  frequency  and  combination  of  modes.  Propaga¬ 
tion  loss  may  be  calculated  by  equations  (32)  and  (33).  The  ranges  over  which  loss 
is  plotted  and  the  source  and  receiver  depths  are  inputted  by  card  group  7.  Values 
of  Dm  are  inputted  by  card  group  8,  and  the  modes  that  make  up  the  sound  field 
are  inputted  by  card  group  9.  A  sample  calcomp  plot  is  shown  in  figure  13. 

The  program  was  evaluated  in  two  ways.  First,  a  test  case  was  run  and  the  re¬ 
sults  compared  to  corresponding  calculations  obtained  from  two  equivalent  procedures 
described  in  references  16  and  17.  Second,  the  convergence  of  the  numerical  pro¬ 
cedure  used  in  solving  the  wave  equation  was  examined  for  many  cases,  and  the 
resulting  limitations  on  the  program  were  noted.  The  test  case  used  for  comparison 
is  illustrated  in  figure  14,  which  gives  the  velocity  profile  and  the  assumed  sound 
velocity  in  the  homogeneous  semi-infinite  bottom. 

Since  each  mode  is  defined  by  a  unique  horizontal  wave  number  Km,  the 
values  of  *cm  obtained  by  each  procedure  were  compared.  Such  a  comparison  is 
shown  in  table  2.  In  the  table,  the  values  of  Km  obtained  by  Bartberger  and 
Ackler  (reference  17)  are  compared  with  the  values  obtained  by  using  the  method 
described  in  this  report.  It  can  be  seen  that  the  results  agree  to  the  fourth  decimal 
place  and,  hence,  tend  to  support  the  validity  of  the  two  procedures.  Furthermore, 
similar  good  agreement  was  obtained  for  this  test  case  using  the  Fast  Field  Program 
(reference  16). 
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Table  2.  Comparison  of  Calculated  *m 


Mode 

xm  (Kanabis) 

*m  (Bartberger  and  Ackler) 

i 

3.78104 

3.78105 

2 

3.77896 

3.77896 

3 

3.77725 

3.77725 

4 

3.77574 

3.77575 

5 

3.77436 

3.77437 

6 

3.77309 

3.77309 

7 

3.77186 

3.77187 

8 

3.77063 

3.77065 

9 

3.76932 

3.76933 

10 

3.76782 

3.76784 

11 

3.76617 

3.76618 

Another  theoretical  result  tested  was  the  calculation  of  propagation  loss  as  a 
functionof  range  for  the  test  case.  The  comparison  of  results  with  those  of  DiNapoli 
(reference  15)  is  shown  in  figures  15  and  16.  It  can  be  seen  from  the  figures  that 
the  interference  patterns  are  virtually  identical  at  all  ranges  except  those  very  near 
the  source.  The  differences  near  the  source  result  from  the  inclusion  by  DiNapoli 
of  "attenuated11  modes  in  the  sound  field.  These  attenuated  modes  contribute  to 
the  field  significantly  at  short  ranges  but  are  attenuated  rapidly  so  that  they  are 
usually  insignificant  at  longer  ranges.  The  procedure  described  in  this  report  con¬ 
siders  only  "unattenuated"  modes,  i.e.,  modes  whose  angles  of  incidence  exceed 
the  critical  angle  so  that  there  is  "total"  reflection  at  the  boundaries.  The  only 
major  discrepancy  in  the  results  lies  in  the  absolute  levels  of  loss  calculated.  For 
this  example,  the  results  of  our  program  consistently  show  3  dB  less  loss  than  the 
values  obtained  from  the  two  other  programs  used  in  the  comparison.  The  source  of 
the  difference  probably  lies  in  the  establishment  of  a  unit  source  level  at  unit  dis¬ 
tance  in  our  procedure. 
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It  was  found  that  the  numerical  procedure  converges  in  practically  all  cases 
considered.  However,  at  high  frequencies,  and  for  extremely  irregular  velocity 
profiles,  the  procedure  may  encounter  difficulty  in  converging.  This  problem  can 
probably  be  overcome  in  most  cases  at  the  expense  of  increased  computer  time, 
e.g.,  by  narrowing  the  interval  between  grid  points  over  the  water  column  (by 
increasing  N  given  in  card  group  5,  table  1). 

It  is  difficult  to  estimate  meaningful  computer  execution  times  because  these 
times  are  sensitive  to  a  large  number  of  input  parameters.  However,  the  order  of 
magnitude  of  time  for  a  typical  case  involving  any  one  of  the  three  options  and  a 
moderate  number  of  incremental  parameters  is  one  minute  on  the  Univac  1108. 


22 


WATER  DEPTH  (ft  ) 

2  00-0  1 80,0  160.0  H  0.0  120  0  100.0  60  0  60.0  40.0  20.0 


TR  4319 


•1.00  -0.50  0-00  0.50  1.00 

AMPLITUDE 


Figure  1.  Amplitude  versus  Depth, 
Frequency  282  Hz,  Mode  1 


Figure  2.  Ray  Equivalent, 
Frequency  282  Hz,  Mode  1 


23 


Figure  3.  Path  of  Reflected  Waves  That  Interfere 
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Figure  4.  Downgoing  and  Upgoing  Waves 
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MODE  I  MODE  1 

PRESSURE  PRESSURE 

- 


Figure  5.  Pressure  of  Upgoing  and  Downgoing  Rays  and  the  Resulfanf 
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BOTTOM 
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Figure  7.  Velocity  Profile,  Complex  Variation 
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Figure  8.  Ray  Equivalent,  Frequency 
200  Hz,  Mode  1 


28 


WATER  DEPTH  (ft) 


TR  4319 


Z--0 


Z 

f 


r 


p .  c , 


Z  =  H 


/>«  c, 

Figure  9.  Two-Layered  Half-Space 
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Figure  10.  Expanded  Velocity  Profile 
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Figure  11.  Sound  Velocity  and  Excitation  Pressure  versus  Frequency,  Mode  1 


Figure  12.  Angle  of  Incidence  versus  Frequency,  Mode  1 
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Figure  13. 


Sample  Calcomp  Plot- 


31 


TR  4319 


BOTTOM 


L 


<900.0  <850  0  5000.0  5O5OJ0  5100.0  51504) 

SOUNO  VELOCITY  (ft/sec) 


Figure  14.  Velocity  Profile,  Test  Case 
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Figure  16.  Propagation  Loss  versus  Range,  Test  Case,  FFP  Program 
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Section  2 

A  COMPUTER  PROGRAM  (S1548)  TO  CALCULATE  NORMAL  MODE 
PROPAGATION  IN  A  MEDIUM  IN  WHICH  STRATIFICATION 
IS  A  FUNCTION  OF  POSITION 


THEORY 

The  normal  mode  theory  used  in  calculating  the  sound  field  in  the  case  of  a 
flat,  homogeneous  ocean  bottom  has  been  discussed  in  section  1  .  The  extension  of 
this  theory  to  a  homogeneous,  semi -infinite  bottom  whose  properties  vary  with  range 
(reference  12)  is  discussed  here. 


In  equation  (34)  of  section  1,  the  pressure  amplitude  of  the  sound  field  produced 
by  a  projector  with  unit  source  level  in  an  ocean  with  a  flat,  homogeneous  bottom 
was  described  by  a  function  dependent  on  range.  The  equation  is  repeated  below. 


Pa  = 


Po 


[e 

|_m 


_n  r/20-ar/20 
PJO  ",  ,  cos(«mr 


Lm 


.  .  -D  r/20-ar/20 

+  I  iL  P^IO  m  sin  ( *mr -tt/4) 


n2 


1/2 


/  (34) 


whe 


re 


a  is  the  attenuation  coefficient  due  to  absorption  in  the  water, 

Po  is  the  water  density, 
r  is  the  horizontal  range, 
m  is  the  mode  number, 

k  is  the  horizontal  wave  number  for  mode  m,  and 
Dm  is  the  attenuation  for  unit  range  due  to  losses  at  the  boundaries. 

Pm  =  Pm  V>m(zo)  ,  (35) 
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wh 


ere 


PQ  (2tt) 


1/2 


pm  2v  ^nr~ 

m  V  m 

+  00 


m 


=  f  P <p  (z)  dz  , 


ps  is  the  water  density  at  the  source,  and 
<Pm(z)  is  the  normalized  displacement  potential  as  a  function  of  depth. 

In  computing  the  sound  field  fora  horizontally  stratified  medium,  the  following 
two  assumptions  (reference  12)  are  made: 

a.  The  values  of  v?m(z)  correspond  to  local  stratification. 

b.  The  stratification  varies  slowly  from  one  region  to  another  so  that  there  is 
no  appreciable  scattering  of  energy  from  one  mode  to  another  when  sound  propa¬ 
gates  through  the  medium. 

If  the  latter  assumption  holds,  then  the  bottom  topography  can  be  approximated 
by  a  number  of  segments  of  different  depth  parallel  to  the  surface  as  shown  in 
figure  17  (reference  18). 


The  power  in  a  given  mode  is  determined  by  <p(zQ),  where  zQ  is  the  source 
depth.  The  pressure  level  at  the  receiver  is  determined  by  <p( z),  where  z  is  the 
receiver  depth.  When  there  is  little  scattering  between  modes,  we  can  approximate 
(reference  19)  the  pressure  level  for  a  given  mode  by  substituting  for  P_  in 
equation  (34),  with 


m 


m 


Pm  -  [P 


P  .  l'/2 

source  receiverj 


(36) 


where  "source"  and  "receiver"  refer  to  the  values  of  P  obtained  by  taking  into 
account  the  local  stratification  at  the  source  or  receiver. 


Thus,  Psource  's  f°ur|d  by  calculating  Pm  in  equation  (35),  using  the  depth, 
velocity  profile,  and  bottom  composition  present  in  the  segment  of  the  range  con¬ 
taining  thesource.  In  addition,  <£m(z0)  should  be  substi  tuted  for  ^  (z) .  Similarly, 
to  calculate  P receiver /  tbe  local  stratification  at  the  receiver  should  be  used  in 
equation  (35)  and  ^  (z)  should  be  substituted  for  <pm (zQ) . 
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In  equation  (34)  the  phase  of  the  signal  depends  upon  *mr.  In  a  horizontally 
stratified  medium,  is  a  function  of  range  since  the  variation  in  water  depth 

and  bottom  composition  causes  to  vary.  The  term  xm  can  be  expressed  by 

(reference  12)  -  _  IfC* 


“mW  =  +  em(r)  / 


r-  > 

kff\ 


m  m ' 


where  is  the  average  of  Km  over  r  and  xm(r)  is  the  value  of  Km  at  a 

range  r.  Equation  (36)  can  be  used  to  form  a  perturbation  solution  of  the  plane 
wave  equation 

2 ,  ^  u2  ^  = 

V  4?  +  — -  <t>  =  0  ; 

c 

in  this  case,  one  obtains  (reference  12)  an  approximate  solution  so  that  in  equation 
(34),  k  (r)  can  be  replaced  by 


Kmr  +  AS 

m  m 


where 


ASm  = 


L 


em(r)  dr  . 


(37) 


If  the  bottom  properties  vary  with  range  from  the  source,  then  the  attenuation 
per  unit  range  due  to  losses  at  the  boundaries  Dm  varies  as  a  function  of  range, 
and  Dm(r)  in  equation  (34)  can  be  replaced  by 


r 

J. 


D  (r)  dr  . 
m 


Thus  equation  (34)  becomes 


Pa 


VT 


„  2 


E  p  io 

m  m 


i 

J  Dm(r)dr/20-ar/20 


cos(  *mr  +  AS  -  it/ 4) 


+ 


Lp  io 

m  m 


•/20-ar/20 


sin(  TT  r  +  AS  -  tt/4) 

m 


J/2 


(38) 
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and 


Lr  =  20  log  Pa  ,  (39) 

where  l_r  Is  the  propagation  loss  at  range  r. 

The  amplitude  distribution  as  a  functlonof  depth  and  the  ray  equivalent  can  be 
calculated  for  a  particular  segment  of  the  range  by  a  method  Identical  to  that  In 
section  1,  using  the  parameters  determined  by  the  local  stratification. 


The  above  formulation  corresponds  to  a  simple  physical  picture  of  the  propaga¬ 
tion  of  sound  In  a  medium  whose  stratification  changes  with  range.  Let  us  consider 
two  plane  waves  propagating  in  a  medium  segmented  as  described  above.  Let  the 
first  wave  correspond  to  a  given  order  mode  In  one  segment  and  the  second  wave 
correspond  to  the  same  order  mode  in  an  adjacent  segment.  If  the  change  In  strati¬ 
fication  between  segments  is  small,  the  difference  in  direction  of  the  two  waves  Is 
small,  and  hence  the  waves  are  excited  at  approximately  the  same  level  at  the 
source.  This  level  is  determined  by  the  excitation  pressure  of  the  mode.  However, 
the  field  at  the  receiver  will  correspond  to  the  stratification  of  the  segment  in  which 
the  receiver  is  located,  and  the  plane  wave  that  corresponds  to  the  mode  considered 
will  dominate  the  contribution  by  that  mode  to  the  field  in  the  segment.  Thus,  the 
field  in  the  segments  is  determined  by  the  plane  wave  corresponding  to  the  modes 
that  compose  the  field.  Extending  the  concept  to  several  segments  and  several  plane 
waves,  we  can  similarly  picture  a  given  mode  changing  shape  in  conformance  to  the 
local  stratification  as  it  moves  from  segment  to  segment. 


Let  us  consider,  for  a  specific  case,  the  effect  upon  sound  propagation  of  dif¬ 
ferences  in  eg,  the  velocity  of  sound  in  the  bottom,  found  in  segments  of  a  range. 
This  effect  can  be  seen  in  table  3,  in  which  9],  the  angle  of  incidence  at  the  bot¬ 
tom  for  the  wave  corresponding  to  mode  1  at  127  Hz,  is  given  as  a  function  of  eg 
for  the  velocity  profile  shown  in  figure  18.  Consider  the  case  where  eg  varies 
between  5200  and  6000  ft/sec  over  a  range.  The  decrease  in  0]  is  rather  small, 
0.9°.  Thus  it  appears  that  the  formulation  presented  above  in  equations  (38)  and 
(39)  can  be  used  for  such  a  variation  in  eg,  with  D](r),  the  bottom  loss  for  mode  1 
as  a  function  of  range,  chosen  to  correspond  to  c^  in  each  segment. 
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Table  3.  Angle  of  Incidence  at  Bottom  for  Mode  1,0], 
and  Critical  Angle  0C  as  a  Function  of  Velocity  of 
Sound  in  Bottom  eg,  at  127  Hz 


CB  (ft/sec) 

91 

0C 

5000 

83°  30' 

82°  22' 

5200 

82°  O' 

72°  22' 

5400 

CO 

o 

CO 

o 

66°  36' 

5500 

V 

CN 

o 

CO 

64°  18' 

5600 

81°  18' 

62°15' 

5800 

81°  12' 

58°  42' 

6000 

o 

o 

CO 

55°  41' 

Let  us  consider  the  result  when  the  value  of  eg  varies  from  5000  to  5200ft/sec. 
The  change  in  0]  is  larger  (1.5°)  than  for  the  case  considered  previously.  Further¬ 
more,  the  0C  corresponding  to  a  eg  of  5000  ft/sec  is  greater  than  the  0]  corres¬ 
ponding  to  5200  ft/sec.  In  the  situation  where  eg  increases  with  increasing  range 
from  5000  to  5200  ft/sec,  the  ray  corresponding  to  mode  1  for  the  harder  bottom 
(eg  =  5200  ft/sec)  can  be  highly  attenuated  because  the  angle  of  incidence  is  less 
than  the  critical  angle  as  the  ray  passes  through  the  segment  that  contains  the  softer 
bottom.  If  the  softer  segment  is  long,  there  is  a  high  probability  that  the  ray  con¬ 
sidered  can  not  contribute  significantly  to  the  sound  field  over  the  hard  bottom. 
Thus,  if  the  ray  is  highly  attenuated,  one  would  expect  a  dramatic  drop  in  signal 
level  and  an  unstable  field.  If,  on  the  other  hand,  the  transition  is  from  a  harder 
to  a  softer  bottom,  the  value  of  0]  increases,  and  hence  the  sound  field  should  be 
stronger  over  the  softer  bottom  because  of  the  transition  alone.  (This  subject  and 
pertinent  experimental  results  are  discussed  in  reference  20.) 

A  change  in  stratification  due  to  an  increase  in  water  depth  results  in  an  in¬ 
crease  in  the  0]  associated  with  a  given  mode.  Hence,  the  result  of  such  a  change 
in  stratification  is  similar  to  the  result  of  a  decrease  in  eg.  (The  effect  of  a  change 
in  depth  is  discussed  in  detail  in  reference  21  in  the  context  of  the  consequences  of 
a  change  in  water  depth  due  to  tide.) 
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Finally,  a  change  in  stratification  that  can  be  described  by  a  change  in 
velocity  profile  with  range  has  two  effects.  First,  like  the  changes  in  strati f i  — 
cation  discussed  above,  the  angle  0]  associated  with  a  given  mode  is  changed. 
Second,  the  angle  of  propagation  of  a  given  wave  corresponding  to  a  given  mode 
in  one  segment  is  diverted  toward  the  direction  corresponding  to  the  same  mode  in 
the  following  segment  entered  by  the  wave. 

DESCRIPTION  OF  PROGRAM  S1548 

Program  Si 548  uses  normal  mode  theory  to  predict  acoustic  propagation  in  a 
medium  over  a  bottom  whose  depth  and  acoustic  properties  vary  with  range.  Such 
a  bottom  is  shown  in  figure  17/ where  the  BIFI  range  is  approximated  by  73  seg¬ 
ments,  each  of  which  has  a  constant  depth.  The  range  is  divided  into  four  regions, 
each  possessing  a  different  bottom  composition.  In  S1548,  provision  is  made  also 
for  the  variation  in  sound  velocity  profile  over  the  range. 

This  section  describes  the  following  aspects  of  the  program: 

a.  Input  parameters  specifying  the  medium  through  which  sound  is  transmitted 

b.  Amplitude  versus  depth  and  the  ray  equivalent 

c.  Propagation  loss  versus  range. 

INPUT  PARAMETERS  SPECIFYING  THE  MEDIUM 

The  number  of  segments  into  which  a  given  acoustic  range  is  divided  is  given 
by  NPCS  (card  group  2,  table  4).  The  depth,  location,  size,  and  bottom  properties 
of  each  segment  are  inputs  of  card  group  8.  The  number  of  velocity  profiles  necessary 
for  a  full  description  of  the  range  is  given  by  NVCP  (card  group  2).  The  number  of 
the  first  segment  described  by  a  particular  velocity  profile  is  given  by  ICVP  (card 
group  3).  This  profile  describes  the  velocity  conditions  from  segment  ICVP(J)  to 
segment  ICVP(J  +  1).  IVPL  (card  group  6)  is  the  parameter  governing  the  option  of 
a  calcomp  plot  of  the  profile. 

AMPLITUDE  DISTRIBUTION  AND  RAY  EQUIVALENT 

Each  segment  of  a  range  has  its  own  amplitude  distribution  as  a  function  of 
depth  for  each  mode.  With  the  range  divided  into  73  segments,  it  is  doubtful  that 
one  would  want  the  amplitude  distribution  and  ray  equivalent  for  every  segment. 

^Figures  17  through  21  oppear  ot  the  end  of  this  section,  poges  45  through  48. 
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Table  4.  Input  Data  for  Program  S1548 


Card 

Amount  of 

Input 

Parameter 

Columns 

Group 

Cards 

Format 

Data 

i 

i 

20A4 

1-80 

Heading  on  each  page  of  output. 

2 

i 

5110 

NPCS 

1-10 

Number  of  segments  of  constant 

depth  into  which  range  is 
divided. 

NVPC 

11-20 

Number  of  velocity  profiles  in 

the  range. 

NADS 

21-30 

Number  of  segments  for  which 
amplitude  distributions  and 
ray  equivalents  are  plotted. 

NCOR 

31-40 

Segment  number  considered 
to  be  the  effective  water 
depth  at  the  source. 

) 

f 

NMOD 

41-50 

Number  of  modes  calculated. 

3 

The 

8110 

ICVP(J) 

1-80 

ICVP(J)  is  the  number  of  the 

smal  lest 

bottom  segment  at  which 

integer 

velocity  profile  J  is  first  used. 

>NVPC 

This  profile  is  used  in  calcu¬ 

8 

lations  between  segments 

ICVP(J)  and  ICVP(J  +  1). 

4 

The 

2110 

IADS(J) 

1-80 

IADS(J)  is  the  number  of  the 

smal  lest 

bottom  segment  for  which  the 

integer 
>  NADS 

pressure  amplitude  distribution 
and  ray  equivalent  are  plotted. 

8 

5 

1 

2110 

IRIC 

1-10 

Range  increments  (ft)  in  propa¬ 

1 

r 

gation  loss  versus  range  plot. 

T 

NPS 

11-20 

Number  of  propagation  loss 
versus  range  plots. 

4F10.3 

ZS 

21-30 

Source  depth  (ft) . 

ZRC 

31-40 

Receiver  depth  (ft). 

FMI 

41-50 

Increment  in  range  in  nmi/in. 
on  propagation  loss  curves. 

* 

FSC 

51-60 

200  x  FSC  is  maximum  depth 
(ft)  plotted  on  curves  that  have 
depth  as  a  parameter. 
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Table  4  (Cont'd).  Input  Data  for  Program  Si 548 


Card 

Group 

Amount  of 
Cards 

Format 

Input 

Parameter 

Columns 

Data 

5 

2F10.3 

FMJ 

61-70 

Value  of  propagation  loss  (dB) 

at  origin  of  propagation  loss 

curves. 

1 

r 

FMK 

71-80 

Increment  in  propagation  loss 

in  dB/in.  on  propagation  loss 

curves. 

6 

Card  a 

110 

NUMV 

1-10 

Number  of  values  in  velocity 

=  i 

profi  le. 

F10.0 

VEL  1 

11-20 

Velocity  at  origin  of  plot  of 

velocity  profile. 

2110 

IVPL 

21-30 

If  IVPL  =  1,  velocity  profile 

1 

not  plotted. 

T 

IEX 

31-40 

Changes  increment  of  kf  by  a 
factor  of  10-IEX.  Values  from 

0  to  10.  (Explained  in  section  1 .) 

Card  b 

2F10.3 

Z  MM 

1-10 

Water  depth  (ft)  at  location  of 

=  1 

j 

s- 

CB 

11-20 

velocity  profile. 

Velocity  (ft/sec)  of  sound  in 

bottom  at  location  of  velocity 
profi  le. 

Card(s)  c 

2F10.3 

ZZ(I) 

1-10 

Height  (ft)  above  bottom  at 

=  NUMV 

\ 

CC(I) 

1 1-20 

which  sound  velocity  is  CC(I). 
Velocity  (ft/sec)  of  sound  at 

ZZ(I). 

1=1,2,. 

NUMV 

in  order  of  increasing  Z. 

(Other 

groups  of  cai 

rds  represe 

1  1 
;nting  a  velocity  profif 

e  should  be  inserted  in  front  of 

Group  8  cards  that  represent  the  bottom  segments  ICVP(J),  J  =  1,  NVPC.) 

7 

1 

110 

N 

1-10 

Number  of  intervals  into  which 

depth  is  to  be  subdivided  in 
integration  of  differential 
equations  (see  section  1). 

2F10.3 

UM 

11-20 

Maximum  value  of  u(z)  (see 

1 

section  1) . 

1 

r 

FQ 

21-30 

Frequency  (Hz), 

42 


TR  4319 


Table  4  (Cont-'d).  Input  Data  for  Program  S1548 


Card 

Group 

Amount  of 
Cards 

Format 

Input 

Parameter 

Columns 

Data 

8 

Card  a 

F10.0 

ZM 

1-10 

Depth  (ft)  of  water  in  a  given 

=  1 

bottom  segment. 

2110 

IRST 

11-20 

Distance  (ft)  of  a  given  bottom 

segment  from  the  source. 

T 

r 

IREN 

21-30 

Length  (ft)  of  a  given  bottom 

segment. 

Card  b 

3F10.3 

CB 

1-10 

Velocity  (ft/sec)  of  sound  in  a 

1 

=  1 

given  bottom  segment. 

RO 

11-20 

Density  (grams/ cm3)  of  water 

above  a  given  bottom  segment. 

1 

1 

RB 

21-30 

Density  (grams/cm3)  of  bottom 

in  a  given  bottom  segment. 

Card(s)  c 

4F10.5 

DD(J) 

1-40 

DD(J)  is  the  loss  at  the  bound¬ 

=  the 

aries  per  unit  range  for  a 

smallest 

given  bottom  segment  (dB/ft 

integer 

for  mode  J) . 

>  NMOD 

J  -  1 1  2,  . .  . ,  NMOD,  where 

4 

NMCD  is  the  number  of  modes 

analyzed . 

(There 

are  NPCS  sets  of  cards 

a,  b,  and  c 

,  where  NPCS  is  the  number  of  segments 

of  constant  depth  into  which  the  range  is  divided.  Set  i  describes  segment  i,  where 

i  =  1, 

21,  NPCS.) 

9 

Card  a 

110 

NMS 

1-10 

Number  of  modes  to  be  summed 

=  1 

in  plot  of  propagation  versus 

range.  NMS  <  NMOD, 

i 

Card(s  b 

6110 

MDS(J) 

1-60 

J  =  1,  2,  .../NMS.  The 

=  the 

values  of  MDS(J)  are  the 

smal  lest 

mode  numi  ers  of  the  modes 

integer 

to  be  summed  in  the  propaga¬ 

>  NMS 

tion  loss  versus  range  plots. 

6 

MDS(J)  <  NMOD. 

(I  here 

are  NPS  'ets  of  Group  9.) 
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Therefore,  NADS  (card  group  2)  is  the  number  of  segments  for  which  calcomp  plots 
of  these  functions  are  desired,  and  IADS(J)  (card  group  4)  gives  the  segment  numbers 
for  which  these  plots  are  generated.  NMOD  (card  group  2)  is  the  number  of  modes 
calculated.  Sample  plots  are  shown  in  figures  19  and  20. 

PROPAGATION  LOSS  CURVES 

Propagation  loss  for  a  source  level  (dB//l  yd)  is  obtained  as  a  function  of  range 
for  any  given  frequency  and  combination  of  modes.  Propagation  lass  may  be  calcu¬ 
lated  from  equations  (38)  and  (39).  Sourceand  receiver  depths  and  ranges  over  which 
loss  is  ta  be  plotted  are  card  group  5  inputs.  Input  of  the  modes  that  make  up  the 
sound  field  is  made  with  card  group  9.  A  sample  plot  is  shown  in  figure  21 . 

It  can  be  seen  from  figure  17  that  the  choice  of  the  effective  water  depth  at 
the  source  may  nat  be  a  straightforward  matter.  The  water  depth  at  the  source  is 
about  17  yd,  ar  51  ft.  However,  most  modes  generated  may  be  influenced  mare  by 
segment  12,  about  37  yd  (110  ft)  deep,  than  by  segments  1  through  11,  since  the 
modes  considered  travel  in  a  direction  characterized  by  large  values  of  the  angle  9 
measured  relative  to  the  normal  to  the  bottom.  If  one  wishes  to  alter  the  effective 
source  depth,  one  would  choose  NCOR  equal  to  the  number  of  the  segment  having 
the  desired  water  depth.  If  one  were  ta  want  the  effective  source  depth  ta  be  that 
af  segment  1,  NCOR  would  be  set  ta  1  .  In  the  case  described,  ane  might  chaose  a 
value  of  NCOR  equal  ta  12.  Propagation  loss  curves  will  be  platted  starting  at 
ranges  corresponding  ta  segment  NCOR. 

No  "unattenuated"  mode  can  possibly  be  formed  at  frequencies  below  the  cut¬ 
off  frequency  for  a  given  mode.  If  this  si  tuati  an  occurs  in  any  segment  of  therange, 
the  propagation  lass  curve  will  be  discontinuous  at  those  ranges  corresponding  ta  the 
location  af  the  particular  segment.  The  effect  af  the  boundaries  represented  by 
/a  Dm(r)dr  aver  this  segment  is  calculated  and  added  to  the  integral  in  equation  (38) 
only  when  it  is  possible  for  "unattenuated"  modes  to  propagate  over  a  range  beyand 
the  range  of  discontinuity. 

SUMMARY 

Program  S1548  is  designed  to  calculate  and  plat  many  quantities  af  interest  in 
the  study  af  a  sound  field  in  an  acean  whose  boundaries  vary  in  depth  and  acoustic 
properties.  Changes  in  the  velocity  profile  aver  the  range  may  also  be  included  in 
the  calculations. 

The  limitations  an  the  velocity  profile  pertinent  ta  program  S1441  (section  1) 
must  be  applied  also  to  program  S1548. 
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Figure  17.  Model  of  BIFI  Range 
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Figure  19.  Amplitude  versus  Depth, 
Frequency  46  Hz,  Mode  1 


Figure  20.  Ray  Equivalent, 
Frequency  46  Hz,  Mode  1 
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Section  3 

PROPAGATION  OF  EXPLOSIVE  SOUND  IN  THE  BIFI  RANGE 

PROCEDURE 

For  all  the  tests  discussed,  explosives  were  detonated  at  depths  of  50  and  75  ft 
at  point  A  in  figure  22/  near  Block  Island.  Signals  were  received  by  a  bottom- 
mounted  hydrophone  located  at  a  depth  of  155  ft  at  point  B,  off  Fishers  Island.  The 
explosive  charges  used  were  1/ 2—  and  1-lb  blocks  of  TNT. 

THEORY 

Most  of  the  theory  of  sound  propagation  by  normal  modes  has  been  discussed  in 
sections  1  and  2.  Section  1  contains  a  description  of  NUSC  program  S1441,  which 
deals  with  normal  mode  propagation  over  a  flat,  homogeneous  bottom  in  a  medium 
whose  velocity  profile  is  constant  with  distance  from  an  acoustic  source.  Section  2 
describes  NUSC  program  S 1 548,  which  uses  normal  mode  theory  to  predict  acoustic 
propagation  over  an  ocean  bottom  whose  depth  and  acoustic  impedance  vary  slowly 
with  range  and  in  a  medium  whose  velocity  profile  varies  slowly  with  distance  from 
an  acoustic  source.  These  two  procedures  will  be  referred  to  as  normal  mode  pre¬ 
dictions  for  a  flat  bottom  and  an  irregular  bottom,  respectively. 

The  prediction  of  propagation  loss  for  a  flat  bottom  is  a  function  of  the  depth, 
velocity  profile,  and  bottom  characteristics,  each  of  which  is  assumed  constant  with 
distance  from  an  acoustic  source.  The  prediction  of  propagation  loss  for  an  irregular 
bottom  is  dependent  on  the  values  of  these  three  parameters  at  both  the  source  and 
receiver .  Thus,  if  the  parameters  have  large  variations  with  distance  from  an  acous¬ 
tic  source,  the  two  methods  may  predict  significantly  di  fferent  values  of  propagation 
loss . 


One  special  case  where  the  predicted  values  do,  indeed,  vary  significantly 
should  be  mentioned.  It  occurs  when  velocity  profi  le  and  bottom  characteristics 
remain  relatively  constant  over  an  entire  range  but  a  large  variation  exists  between 
the  depth  assumed  for  a  flat  bottom  and  the  depth  at  the  source  or  receiver.  Vertexing 
of  the  rays  can  account  for  a  marked  disparity  between  the  sound  fields  at  the  re¬ 
ceiver.  As  shown  in  figure  23,  H  ^  is  the  assumed  depth  for  a  flat  bottom  case  and 
H2  is  the  depth  at  the  source  or  receiver  in  the  irregular  bottom  case;  the  source  or 
receiver  is  at  the  bottom.  In  the  irregular,  but  not  the  flat,  bottom  case  the  hydro¬ 
phone  is  in  a  "shadow  zone"  and  the  calculated  propagation  loss  is  consequently 
* 

Figures  22  through  37  appear  at  the  end  af  this  section,  pages  60  through  73. 
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much  larger  since  the  omplitude  of  the  pressure  fieldsshows  on  exponential  decoy 
in  o  shadow  zone  (reference  8).  The  example  considered  occurs  in  the  winter  ond 
where  H  -j  =  1 1 0  ft  is  the  overage  depth  used  for  o  flat  bottom  ond  H2  =  150  ft  is 
the  depth  at  the  receiver. 

It  has  been  previously  determined  (references  7  ond  8)  thot  ot  the  frequencies 
considered  in  these  tests  (about  56  to  560  Hz),  the  first  mode  dominates  in  the 
signals  received  ot  point  B  in  figure  22.  Therefore,  it  hos  been  ossumed  in  the 
theoretical  colculotions  of  propagation  loss  thot  the  pressure  field  at  the  receiving 
hydrophone  contains  only  the  first  mode.  Since  the  first  mode  is  dominant,  there  is 
little  difference  in  the  pressure  field  produced  by  explosives  ond  detected  ot  50  ond 
75  ft.  As  in  references  7  ond  8,  the  measurements  ot  the  two  depths  ore  combined. 
Theoretical  calculations  ore  for  o  source  75  ft  deep. 

Let  us  consider  the  general  chorocteristics  of  propagation  loss  os  o  function  of 
frequencyin  cases  where  the  first  mode  dominates.  All  calculotions  are  made  assuming 
the  velocity  profile  shown  in  figure  24.  Figure  1 1  is  o  typicol  plot  of  excitation  pres¬ 
sure  (defined  in  section  1  of  this  report)  versus  frequency  for  the  first  mode.  At  the  fre¬ 
quencies  considered  here,  excitation  pressure  decreases  with  frequency.  Thus,  less 
energy  goes  into  the  higher  frequency  modes  and,  because  of  this  factor  olone,  propa¬ 
gation  loss  wi  1 1  increase  with  increasing  frequency.  It  is  possible  (reference  8)tocon- 
struct  a  ray  equivalent  of  a  particular  mode,  os  shown  in  figure  25;  it  is  therefore 
possible  to  determine  the  skip  distance  between  bounces  off  the  bottomof  the  medium. 
The  angle  at  which  energy  strikes  the  surface  0$  or  bottom  0^ ,  relotive  to  the 
normal  to  the  bottom,  increases  with  frequency  (reference  8) .  Thus,  as  shown  in  figures 
25  ond  26,  the  skip  distance  tends  to  increase  with  frequency,  thereby  reducing  the 
number  of  bounces  over  o  given  ronge.  However,  if  a  velocity grodient  exists  in  the 
medium,  the  skip  distance  increases  with  the  frequency  until  there  is  a  vertexing  ot 
the  interface,  as  presented  in  figure  27.  In  generol,  this  represents  the  lorgest  skip 
distance  attainable .  As  the  frequency  is  increased  further,  the  depth  of  the  vertexing 
recedes  from  the  interface  ond  the  skip  distonce  decreases,  os  shown  in  figure  28. 
Thus,  propagation  loss  wi  1 1  decrease  wi  th  increasing  frequency  as  o  result  solely  of  the 
effect  of  skip  distonce;  however,  the  effect  of  vertexing  con  modify  this  relationship. 

In  the  irregulor  bottom  cose,  mony  depths  and  velocity  profiles  moy  exist. 
Obviously,  then,  there  are  many  possible  ways  todetermine  the  skip  distance.  One 
woy  is  totoke  on  averoge  skip  distonce  for  thevorious  sections  of  an  acoustic  ronge. 
However,  the  sections  of  the  range  with  the  largest  gradients  should  produce  the 
greotest  effect  on  the  propagation  loss.  Therefore,  in  the  analysis  of  these  tests, 
whenever  the  velocity  profile  varied  considerably  over  the  ronge,  the  skip  distonce 
was  determined  by  using  the  largest  grodient  measured  ot  o  depth  thot  was  close  to 
the  averoge  depth  of  the  range. 
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RESULTS  OF  TESTS 

Propagation  loss  as  a  function  of  frequency  was  determined  for  all  five  tests 
conducted.  The  results  were  obtained  by  finding  the  energy  content  of  each  re¬ 
ceived  shot  for  a  1-Hz  band  at  logit  frequencies  from  56  to  562  Hz.  The  levels  thus 
derived  were  subtracted  from  the  source  level  for  the  explosives  as  given  by  Weston 
(reference  22). 

For  each  test,  the  theoretical  propagation  loss  values  were  calculated  in  both 
the  irregular  and  flat  bottom  cases  with  the  assumption  that  no  bottom  loss  was  suf¬ 
fered  by  the  first  mode .  Velocity  profi  les  measured  at  the  time  of  the  tests  were  used  . 
The  profiles  obtained  during  the  five  tests  ranged  from  those  with  large  negative 
gradients  to  those  with  small  positive  gradients.  In  the  flat  bottom  case,  the  veloc¬ 
ity  profile  with  the  largest  gradient  was  chosen  to  apply  to  the  whole  range  if  more 
than  one  profile  was  taken .  The  water  depth  assumed  for  the  flat  bottom  was  110  ft. 
The  difference  between  experimental  and  theoretical  values  of  propagation  loss  was 
interpreted  as  a  measure  of  bottom  loss,  and  the  internal  consistency  of  the  theoret¬ 
ical  and  experimental  results  was  observed .  The  outcomes  of  the  individual  tests  are 
given  below. 

AUGUST  1967 

The  August  1967  tests  were  conducted  when  the  velocity  profile  possessed  a 
large  negative  gradient,  as  shown  in  figure  24.  Since  this  was  the  only  velocity 
profile  taken  during  these  tests,  it  was  used  to  represent  velocity  conditions  over 
the  entire  range.  As  shown  in  figure  29,  propagation  loss  as  a  function  of  frequency 
is  nonlinear,  with  the  minimum  at  141  Hz.  The  results  and  the  theoretical  analysis 
aregiven  in  table  5.  The  differences  between  theoretical  predictions  of  propagation 
loss  for  the  flat  and  irregular  bottom  cases  are  small,  since  in  both  cases  only  one 
velocity  profile  was  used  to  represent  the  entire  range.  The  increase  in  theoretical 
loss  with  frequency  is  a  result  mainly  of  the  variation  of  the  excitation  pressure 
with  frequency.  At  low  frequencies,  the  skip  distance  increases  with  frequency,  as 
shown  in  figures  25  and  26.  This  increase  continues  until  a  frequency  is  reached  at 
which  the  ray  vertexes  near  the  surface  (figure  27).  This  event  corresponds  to 
141  Hz  in  table  5,  at  which  point  the  skip  distance  is  maximum  and  the  angle  at 
which  the  energy  first  strikes  the  surface  becomes  90°.  Thereafter,  as  shown  in 
figure  28,  the  ray  vertexes  at  increasing  depth  for  higher  frequencies  and  the  skip 
distance  decreases.  The  loss  per  nautical  mile  is  minimum  at  141  Hz,  and  the  loss 
per  bounce  is  on  the  order  of  0.5  dB  over  the  frequency  range  considered.  The 
maximum  loss  per  bounce  was  determined  to  occur  at  141  Hz,  which  is  probably 
a  result  of  the  assumption  that  the  particular  velocity  profile  measured  is  a  good 
approximation  at  all  points  along  the  range. 
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Table  5.  Results  of  August  1967  Tests 
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JANUARY  1968 


The  tests  in  January  1968  were  conducted  when  the  typical  velocity  profile 
had  a  small  positive  gradient.  All  five  profiles  taken  over  the  range  were  similar 
to  that  shown  in  figure  30  (taken  from  reference  23).  For  comparison,  figure  31  ex¬ 
hibits  pi  ots  of  propagation  loss  as  a  functionof  frequency  for  both  the  January  1968 
and  the  August  1967  tests.  It  can  be  seen  that  propagation  loss  was  much  greater  in 
August  than  in  January,  and  that  the  difference  in  propagation  loss  between  the 
two  sets  of  tests  is  particularly  large  for  frequencies  above  200  Hz.  The  results  and 
theoretical  analysis  pertaining  to  the  January  tests  are  shown  in  table  6.  It  can  be 
seen  that  the  predicted  propagation  loss  values  for  an  irregular  bottom  are  much 
larger  than  those  for  a  flat  bottom,  especially  at  the  higher  frequencies.  This  is 
attributed  to  the  previously  described  effect  caused  by  irregularity  in  depth  near 
the  receiver  and  a  positive  velocity  gradient.  In  both  cases,  the  predicted 
propagation  loss  increases  with  frequency,  as  expected.  The  skip  distance  in¬ 
creases  with  frequency  since  vertexing  would  not  occur  until  0l  equals  90° 
at  about  700  Hz. 


It  should  be  noted  that  the  physical  pictures  suggested  by  the  analysis  of  a  flat 
and  an  irregular  bottom  differ.  For  the  flat  bottom,  the  loss  per  bounce  is  on  the 
order  of  0.2  dB  at  all  frequencies  considered.  However,  theskip  distance  increases 
(number  of  bounces  decreases)  with  frequency,  and  thus  bottom  loss  decreases  with 
frequency,  as  demonstrated  by  values  of  bottom  loss  in  dB  per  nautical  mile. 


In  the  irregular  bottom  analysis,  the  loss  per  bounce  decreases  sharply  with 
increasing  frequency  as  the  angle  of  incidence  of  energy  striking  the  bottom,  9^, 
increases  from  74.4  to  89.3°.  This,  when  combined  with  the  increase  in  skip  dis¬ 
tance  with  frequency,  explains  the  decreased  bottom  loss  with  frequency.  At  fre¬ 
quencies  of  355  to  562  Hz,  the  predicted  values  of  propagation  loss  are  lower  than 
the  measured  values  by  about  2  dB,  and  bottom  lossand  loss  per  bounce  are  negative. 
At  low  frequencies,  theskip  distances  for  August  and  January  are  about  the  same. 
However,  the  angle  of  incidence  is  larger  during  January,  a  fact  that  explains  the 
decrease  in  propagation  loss  in  January  at  these  frequencies.  At  higher  frequencies, 
both  the  angle  of  incidence  and  the  skip  distance  are  larger  in  January  than  in 
August.  Thus,  the  larger  seasonal  increase  in  propagation  loss  at  the  higher  fre¬ 
quencies  is  explained. 
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Table  6.  Results  of  January  1968  Tests 
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Number 

of 

Bounces 

141 

118 

96 

80 

65 

54 

44 

35 

28 

19 

18 

Irregular  Bottom 

Loss  Per 
Bounce 

(dB) 

.24 

.20 

.17 

.11 

.11 

.12 

.07 

.04 

-.06 

-.11 

-.12 

Bottom 

Loss 

(dB/  nmi) 

1.86 

1.32 

0.91 

0.48 

0.38 

0.35 

0.17 

0.07 

-0.09 

-0.12 

-0.11 

Flat  Bottom 

Loss  Per 
Bounce 
(dB) 

.26 

.23 

.21 

.16 

.17 

.21 

.19 

.20 

.17 

.27 

.39 

Bottom 

Loss 

(dB/nmi) 

2.01 

1.48 

1.09 

0.69 

0.61 

0.60 

0.45 

0.38 

0.26 

0.28 

0.38 

Measured 

Loss 

(dB) 

93.3 

85.3 

80.2 

75.0 

75.9 

78.3 

78.3 

79.9 

80.9 

84.8 

90.6 

Theoretical  Loss 

Irregular 

(dB) 

59.2 

61.2 

63.6 

66.2 

68.9 

71.9 

75.1 

78.6 

82.5 

87.0 

92.7 

Flat 

(dB) 

56.6 

58.3 

60.3 

62.4 

64.7 

67.3 

70.0 

73.0 

76.1 

79.6 

83.6 

,„'a) 

®  3 

74.1 

76.5 

78.8 

80.5 

82.0 

83.3 

84.3 

85.1 

85.8 

86.2 

86.3 

o  O) 

CD  V 
~o 

74.4 

76.7 

79.2 

81.2 

82.8 

84.3 

85.5 

86.5 

87.5 

88.5 

89.3 

Ski  p 
Distance 

(fO 

780 

930 

1140 

1380 

1680 

2050 

2520 

3120 

3980 

5640 

6160 

Frequency 

(Hz) 

56 

70 

89 

112 

141 

178 

224 

282 

355 

446 

562 
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APRIL  1968 

The  tests  in  April  1968  were  conducted  when  the  velocity  profile  varied  con¬ 
siderably  over  the  range,  as  shown  in  figure  32  (taken  from  reference  23),  As  can  be 
seen,  the  profi le  possessed  a  small  negative  gradient  approximately  18  nmi  from  the 
source  and  a  fairly  large  negative  gradient  in  the  middle  of  the  range.  For  compar¬ 
ison,  plots  of  propagation  loss  as  a  function  of  frequency  for  the  April  and  January 
tests  are  shown  in  figure  33.  It  is  apparent  that  the  propagation  loss  at  most  frequen¬ 
cies  was  slightly  greater  in  Apri  I  than  in  January .  The  April  results  and  the  theoretical 
analysis  are  presented  in  table  7.  The  di fferences  between  the  theoretical  propaga¬ 
tion  loss  predictions  for  a  flat  and  an  irregular  bottom  are  small  since  there  is  little 
variation  in  the  velocity  profile  near  the  source  and  receiver.  As  expected,  the 
theoretical  propagation  loss  increases  with  frequency .  The  skip  distances  and  angles 
of  incidence  at  the  surface  and  bottom  were  calculated  for  about  12  nmi  from  the 
source,  using  the  velocity  profile  shown  in  figure  32.  By  comparing  the  values  of 
skip  distance  and  0^  in  tables  6  and  7,  one  sees  that  these  quantities  are  nearly 
the  same  in  January  and  April  for  low  frequencies,  so  that  similar  bottom  losses 
would  be  expected.  However,  at  higher  frequencies,  the  skip  distance  and  0^  are 
greater  in  January,  so  that  one  would  expect  slightly  lower  values  of  propagation 
loss  in  January.  This  seems  to  be  the  case,  although  the  differences  at  112  and 
141  Hz  are  larger  than  anticipated.  The  loss  per  bounce  in  the  higher  frequency 
range  was  on  the  order  of  0.3  dB. 

AUGUST  1968 

The  tests  in  August  1 968  were  conducted  when  the  velocity  profi les*  taken  over 
the  range  exhibited  large  negative  gradients,  as  shown  in  figure  34  (taken  from  refer¬ 
ence  23) .  The  profi  le  near  the  source,  however,  possessed  only  a  slightly  negative  gra¬ 
dient,  smaller  thon  thegradient  in  the  profi  le  observed  during  the  August  1967  tests. 

Figure  35  gives  plots  of  propagation  loss  as  a  function  of  frequency  for  the 
August  1968  and  the  August  1967  tests.  It  can  be  seen  that  the  propagation  loss  at 
most  frequencies  was  slightly  greater  in  1968  than  in  1967.  These  results  and  the 
theoretical  analysis  for  the  August  1968  tests  are  provided  in  table  8.  There  are 
only  moderate  differences  between  the  theoretical  predictions  of  propagation  loss 
fora  flat  and  an  irregular  bottom,  since  there  is  only  a  moderate  difference  between 
the  velocity  profi  les  at  thesource  and  at  the  receiver .  Again,  the  theoretical  propa¬ 
gation  loss  increases  with  frequency. 

The  skip  distance  is  maximum  at  178  Hz,  compared  with  themaximum  skip  dis¬ 
tance  at  141  Hz  in  August  1967,  In  August  1968,  the  minimum  propagation  loss  at 
141  Hz  is  slightly  less  pronounced  than  the  minimum  loss  for  August  1967. 

*These  profiles  were  token  a  doy  after  the  tests  were  performed. 
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Table  7.  Results  of  Apri I  1968  Tests 
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Table  8.  Result's  of  August  1968  Tests 
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SEPTEMBER  1968 

The  September  1968  tests  were  performed  when  the  velocity  profiles*  taken 
over  the  range  exhibited  moderately  negative  gradients,  as  shown  in  figure  36  (taken 
from  reference  23).  The  profiles  near  the  source  and  receiver  were  less  negative  than 
those  toward  the  center  of  the  range.  Figure  37  gives  plots  of  propagation  loss  as  a 
function  of  frequency  for  the  September  1968  and  August  1967  tests.  The  loss  is  lower 
in  August  at  low  frequencies  and  lower  in  September  at  the  higher  frequencies. 
Table  9  shows  that  in  September  the  maximum  skip  distance  occurs  at  355  Hz;  in 
August  the  skip  distance  is  maximum  at  141  Hz.  Since  the  skip  distances  are  longer 
in  August  at  lower  frequencies  and  longer  in  September  at  higher  frequencies,  it  is 
not  surprising  that  relatively  less  loss  occurs  at  the  higher  frequencies  in  September 
and  at  the  lower  frequencies  in  August. 


*These  profiles  were  abtoined  two  days  after  the  tests  were  performed. 
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Table  9.  Results  of  Sept-ember  1968  Tests 


TR  4319 


L. 

on 

9 

0) 

c 

u-  ^ 

CN 

CN 

no 

CN 

Os 

X 

O 

CO 

N- 

CN 

00 

O  D 

• — 

Os 

x 

O 

N- 

CO 

CO 

CN 

co 

O 

NO 

7 

O 

CQ 

L_ 

E 

CL.  <->  — 

o 

C  CD 

O 

CN 

CO 

O 

CN 

X 

LO 

X 

00 

Os 

-4— 

o 

on  -5  -7-5 
on  ^ 

O  J? 

N- 

N- 

N- 

NO 

Os 

O 

N- 

X 

N- 

N- 

CD 

_l  00 

' 

' 

o 

*"  ' 

D 

O 

|  «.  e 

X 

CN 

00 

O 

, — 

NO 

00 

Os 

0 

CN 

O 

<D 

L_ 

O  on  r- 
£  °  \ 

N- 

• 

x 

LO 

NO 

00 

X 

X 

N- 

O 

00 

1- 

O  — 1  CD 

CN 

CN 

1 — 

1 — 

1 — 

, — 

, — 

1 — 

, — 

, — 

, — 

co  -a 

N — 

!_ 

0)  ^ 

CO 

O 

00 

_ 

O 

X 

0 

0 

CO 

CO 

CN 

E 

o 

s  gs 

0  ° 

N- 

• 

N- 

N- 

10 

X 

0 

CO 

00 

0 

0 

NO 

-4— 

_J  CQ 

' 

o 

CD 

^ _ v 

O 

g  «.  E 

CO 

O 

Os 

00 

00 

XT 

CN 

0 

LO 

O 

LL. 

i  °\ 

0 

CO 

O 

X 

00 

■ 

CN 

Os 

CO 

O  — 1  CD 
CD  "O 

CN 

CN 

1 — 

' — 

CN* 

CN 

CN 

CN* 

CN 

"a 

<u 

j 

— 

on  x — 

LO 

CN 

CN 

00 

CN 

NO 

CN 

LO 

NO 

X 

0 

on  oD 

• 

• 

• 

• 

• 

« 

• 

• 

# 

o 

9  3. 

LO 

- - 

IX 

xf 

00 

LO 

NO* 

00 

CO 

00* 

CN 

<u 

_ 1 s- — 

0 

O 

Os 

CN 

CN 

O 

O 

0 

0 

O 

i/) 

«/> 

i— 

o 

_o 

_i 

2,  ST 

CO 

n- 

>0 

O 

00 

— 

00 

O 

00 

o 

n? 

0 

CN 

N- 

X 

Os 

, — 

N- 

o’ 

X* 

CN* 

CN 

u 

<D  ' 

L- 

0 

0 

N0 

0 

NO 

X 

X 

X 

X 

X 

X 

-4— 

<u 

L_ 

— 

o 

0) 

-4—  /-v 

O  CD 

N- 

1 

00 

0 

N- 

0 

^r 

CO 

00 

• 

CO 

X 

_c 

lZ  3- 

IX 

Os’ 

O 

CN* 

N-* 

0* 

X* 

00* 

00 

Os* 

R 

U-) 

u-> 

O 

NO 

NO 

NO 

0 

NO 

NO 

0 

Os 

CN 

CN 

CX 

NO 

CO 

0 

O 

0 

0 

O 

CD 

x 

O 

CN 

CO 

UD 

X 

0 

O 

0 

0 

O 

x 

00 

00 

00 

00 

00 

Os 

Os 

Os 

CN 

CN 

_Q  ^ 

0 

sO 

00 

00 

NO 

NO 

Os 

CO 

NO 

®  -S 

0 

00* 

Os* 

O 

, — 

CN 

CO 

CO* 

CO* 

3 

N- 

v — 

x 

IN 

IX 

00 

00 

00 

00 

00 

CO 

00 

O 

( 

a. 

U 

C  ^ 

0 

O 

O 

O 

0 

O 

0 

0 

0 

0 

O 

O 

00 

O 

N" 

X 

CO 

ID 

CO 

28 

X 

X 

CN 

-V 

cr) 

S-  o+_ 

Os 

• — 

N" 

X 

CN 

Os 

0 

CN 

X 

NO 

.12 

• — 

1 — 

’ — 

CN 

CN 

CO 

N- 

CO 

1 — 

Q 

X 

u 

! 

Z  ✓ - V 

0)  N 

p 

Os 

CN 

r— 

00 

N- 

CN 

IT) 

$ 

CN 

X 

3  T 

IN 

00 

• — 

N- 

X 

CN 

00 

uo 

NO 

O 

cr  — ^ 
o 

1 

' 

1 

CN 

CN 

CO 

uo 

X 

j 

u_ 

59 


TR  4319 


DISTANCE  FROM  SOURCE  (nmi) 


Figure  22.  Depth  Profile,  BIFI  Range 


Figure  23.  Effect  of  Variable  Depth  on  Rays 
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Figure  24.  Velocity  Profile,  8  August  1967 
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Figure  25.  Ray  Equivalent,  Frequency 
70  Hz,  Mode  1 
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Figure  26.  Ray  Equivalent,  Frequency  112  Hz,  Mode  1 
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Figure  27.  Ray  Equivalent, 
Frequency  141  Hz,  Mode  1 


Figure  28.  Ray  Equivalent, 
Frequency  224  Hz,  Mode  1 


64 


TR  4319 


FREQUENCY  (Hz) 

Figure  29.  Propagation  Loss  versus  Frequency,  August  1967 
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Figure  31.  Propagation  Loss  versus  Frequency,  January  1968  and  August  1967 
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Figure  33.  Propagation  Loss  versus  Frequency,  January  and  April  1968 
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Figure  35.  Propagation  Loss  versus  Frequency,  August  1967  and  August  1968 
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Figure  37.  Propagation  Loss  versus  Frequency,  August  1967and  September  1968 
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Section  4 

SUMMARY  AND  CONCLUSIONS 

Program  S1441  is  designed  to  calculate  and  plot  many  quantities  of  interest  in 
the  study  of  a  sound  field  in  an  ocean  bounded  by  flat  parallel  boundaries.  The 
problem  of  acoustic  propagation  is  approached  from  the  standpoint  of  physical  and 
ray  acoustics. 

Program  S1548  is  designed  to  calculate  and  plot  many  quantities  of  interest  in 
the  study  of  a  sound  field  in  an  ocean  whose  boundaries  vary  in  depth  and  acoustic 
properties.  Changes  in  the  velocity  profile  over  the  range  may  also  be  included  in 
the  calculations. 

The  theoretical  predictions  using  programs  S1441  and  S1548  differ  by  small  or 
moderate  amounts  in  four  of  the  explosive  tests  considered.  In  the  January  1968 
tests,  the  difference  is  considerable,  and  predictions  using  program  S1548  are 
physically  more  plausible .  This  is  true  also  to  a  lesser  degree  in  the  other  explosive 
tests.  Therefore,  assumption  of  a  medium  in  which  stratification  is  a  function  of 
position  seems  to  be  an  improvement  over  that  of  a  medium  wi  th  constant  stratification 
wi th  range . 

The  results  derived  from  the  explosive  tests  are  consistent  wi  th  normal  mode  pre¬ 
dictions  .  Three  major  factors  account  for  the  relationship  between  propagation  loss 
and  frequency.  First,  excitation  pressure  decreases  with  increasing  frequency,  which 
has  the  effect  of  increasing  propagation  loss  with  increasing  frequency.  Second,  skip 
distance,  in  general,  increases  with  frequency,  thereby  decreasing  propagation  loss 
with  increasing  frequency.  In  the  BIFI  range,  these  two  effects  seem  to  cause  the 
minimum  propagation  loss  at  a  frequency  around  100  to  200  Hz.  This  minimum  is 
either  enhanced  or  depressed  by  the  third  factor,  the  frequency  at  which  vertexing 
commences.  If  vertexing  occurs  near  the  minimum,  as  in  August  1967  and  August 
1968,  the  minimum  is  enhanced;  if  it  occurs  away  from  the  minimum,  the  minimum 
is  rendered  less  pronounced. 

An  interesting  effect  explained  by  the  normal  mode  analysis  is  the  dependence 
of  propagation  loss  on  the  size  of  the  negative  gradient  of  the  velocity  profile.  In 
general,  propagation  loss  will  increase  with  thesizeof  the  negative  gradient  because 
an  increase  in  the  negative  gradient  will  tend  to  decrease  the  angle,  0^,  at  which 
energy  strikes  the  bottom,  thereby  decreasing  the  skip  distance.  However,  this  in¬ 
crease  in  negative  gradient  also  lowers  the  frequency  at  which  vertexing  first  occurs, 
which  normally  corresponds  to  the  largest  skip  distance  at  any  frequency.  So,  in 
August  1967  and  August  1968,  propagation  loss  at  the  lower  frequencies  is  less  than 
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the  corresponding  loss  in  September  1968,  even  though  much  larger  negative  gradients 
were  observed  in  the  velocity  profiles  taken  during  the  August  tests.  This  circum¬ 
stance  may  be  attributed  to  the  fact  that  vertexing  took  place  at  about  150  Hz  in 
August  and  at  about  350  Hz  in  September.  The  August  skip  distances  at  low  fre¬ 
quencies  were  larger  than  the  corresponding  skip  distances  in  September.  This  is 
compatible  with  normal  mode  theory,  which  states  that  at  low  frequencies,  propa¬ 
gation  loss  should  increase  to  maximum  as  the  gradient  becomes  negative  and  then 
decrease  with  a  further  increase  in  the  size  of  the  negative  gradient.  At  the  higher 
frequencies  in  the  case  of  the  profiles  considered,  propagation  loss  increases  with 
increasingly  negative  gradients. 


The  determination  of  the  exact  skip  distance  is  a  major  problem  in  normal 
mode  analysis.  As  explained  previously,  the  skip  distance  assumed  over  the  entire 
range  was  calculated  by  using  the  velocity  profile  with  the  largest  gradient.  As  a 
consequence,  whenever  there  is  a  large  variation  in  the  velocity  profile  over  the 
range,  there  is  a  bias  in  the  calculations  of  skip  distance  as  a  functionof  frequency. 
This  can  be  seen  in  table  9,  where  the  calculated  loss  per  bounce  is  extremely  high 
at  355  Hz,  the  frequency  at  which  the  skip  distance  is  maximum.  The  bias  is 
especially  severe  when  the  bottom  loss  is  high  and  the  velocity  profile  varies 
significantly  over  the  range. 

An  a  I  te  motive  method  of  determi  ni  ng  ski  p  distance  wou  Id  be  to  take  an  average 
of  the  skip  distance  over  each  segment  as  weighted  by  the  length  of  each  segment. 
Onedrawback  to  this  method  is  that  segments  with  large  gradients  would  be  weighted 
evenly  with  those  possessing  small  gradients;  the  effect  of  the  smal I—  gradient  seg¬ 
ments  on  propagation  loss  might  be  smaller  than  the  effect  of  the  large-gradient 
segments.  If  an  analysis  similar  to  the  one  described  here  is  performed  in  the  future, 
program  S1548  should  be  modified  to  do  the  calculations  suggested,  and  the  results 
should  be  compared  with  skip  distances  calculated  for  the  segment  with  the  largest 
velocity  gradient.  Ideally,  an  investigator  would  want  to  know  the  bottom  loss  per 
bounce  in  each  segment  so  that  he  could  determine  the  total  loss  by  summing  the 
product  of  loss  per  bounce  by  the  number  of  bounces  in  each  segment.  However, 
experimental  determination  of  bottom  loss  as  a  function  of  frequency  for  angles  of 
about  75  to  90°,  relative  to  the  normal,  would  be  extremely  difficult  to  perform. 


As  expected,  it  was  found  that,  for  a  given  frequency,  loss  per  bounce  de¬ 
creased  as  the  angle  of  incidence  relative  to  the  normal  increased.  Loss  per  bounce 
as  a  function  of  frequency  did  not  show  any  marked  trends  except  for  the  January 
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data  in  the  case  of  an  irregular  bottom.  Here,  the  angle  of  incidence  increased 
extremely  rapidly  with  frequency,  and  the  bottom  loss  went  to  zero  as  the  rays 
nearly  vertexed  close  to  the  bottom.  Since  the  angle  of  incidence  at  the  bottom 
increases  with  frequency,  it  might  be  concluded  in  the  other  cases  that  the  bottom 
loss  for  a  constant  incident  angle  would  also  increase  with  frequency. 

The  theoretical  predictions  for  flat  and  irregular  bottoms  differ  by  small  or 
moderate  amounts  in  four  out  of  the  five  tests  considered.  In  the  January  1968  tests, 
the  difference  is  considerable,  with  the  predictions  for  the  irregular  bottom  physi cal  ly 
more  plausible.  This  is  true  also  to  a  lesser  degree  in  the  other  tests.  Therefore, 
the  assumption  of  an  irregular  bottom  seems  to  be  an  improvement  over  the  assump¬ 
tion  of  a  flat  bottom. 

Further  investigation  of  the  effects  described  in  this  report  will  be  conducted 
during  daily  propagation  tests  at  frequencies  of  127,  400,  and  1702  Hz. 
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Appendix  A 

LIST  OF  PROGRAM  S1441 

DOUBLE  precision  f»u,up 

DOUBLE  PRECISION  K2 ,  W2,P2,H,ZP,  CZ.G,H2,0EN 

DOUBLE  PRECISION  S , K20 ,DK2 , GAM, BA 

DOUBLE  PRECISION  E,V 

DOUBLE  PRECISION  UMAX 

DOUBLE  PRECISION  C0S2  ,  CCC  t  CV 

DIMENSION  hED(20> »Z(lOO) ,C(100) ,F(500) ,UP<500)  ,U<500)  * 

1  ZMP ( 1000 ) •  XMP ( 1000 ) t  C0S2 ( 500 ) ,  CCC(500) 

1  •  ZZP(500)  *  ZPP(500),  VV  (500 )  ,  DATAU024)  »  UU(500) 

I  » PM ( 3000 )  »DDT(100) 

1  »  ZR <500 )  »  CR<500)  »RP<3000> ,ZBT(5>  ,  ZBF(5) 

1 » DU ( 100 )  ,  PPT ( 3000 )  »PCT  <3000)  , PST (3000)  t  MDS ( 100  )  ,USR  < lOo )  »URR(100)  , 
1  PMM <100)  , FMR ( 100 )  »FDS<100) 

1  »GVEL<2000>»  PVEL  <2000 ) ,  FQP<2000),  Tl<2000),  T2<2000) 

COMMON  MS,JM,U 
COMMON  UMS 

CALL  PLOTS  (  UATA ( 1 )  *  1024,  6  ) 

4  JPM  =  25 
READ  <3 , 1 ) HED 

1  FORMAT  (20A4) 

READ  (3,112)  NUMV,  VEL1  t  IVOP  ,IRP  ,  IVPL  ,  lEX 

112  FORMAT  (  110, F10. 0,4110) 

NUM  =  0 

READ (3,2)  ZM,Cd,RO»RB  ,FSC 

2  FORMAT  (5F10.3) 

ZMb  =  (200 , 0*FSC  .  -ZM)/(20.0*FSC) 

ZBT(l)  =  ZM 

ZbT (2)  =  ZM 

ZbT (3)  =  200.0  *FSC 

ZbT (4)  =  -z0.0*FSC 

ZBF(l)  =  0.0 

ZBP (2)  =  10.0 

ZBF (3)  =  0.0 

ZBF  <  4 )  =  1.0 

DO  3  1=1,100 

Z(l)  =  0. 

3  cm  =  o. 

1=1 

5  READ  (3,30)  Z(I) »C(I) 

30  FORMAT (2F10. 3) 

IF (DABS (Z ( I ) -ZM) 01 ) 6 ,6 » 7 
71=1+1 
60  TO  5 

6  1  =  1 

IF  (IVPL. EO.l)  60  TO  8 
DO  130  K  =  1,  NUMV 

CR (NUMV  -K  +1)  =  C ( K ) 

130  ZR (NUMV  -K  +1)  =  Z(K) 

ZR ( NUMV  +1)  =  zM  -0,01 
ZR ( NUMV  +2)  =  ZM  +50.0 
CR ( NUMV  +1)  =  Cd 

CR (NUMV  +2)  =  Cd 

DO  165  N1NY  =  1,  NUMV 
lb5  ZR(NINY)  =  ZM  -  ZR(NXNY) 

ZR (NUMV  +3)  =  200.0+FSC 
ZR (NUMV  +4)  =  -20.0*FSC 
CR ( NUMV  +5)  =  VEL1 
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CRlNUMV  +4)  =  50,0 

CALL  PLOT  (0,0, 0.0, -3) 

CALL  LINE  (CR»ZR»NuMV  +2#  1,0»0  ) 

CALL  LINE  (ZBF,ZBT,2,1,0,0) 

CALL  SYMBOL  (10.25,ZMB,0.14,6HBOTTOM»0.0#6  ) 

CALL  AXIS  (0.0,0,0,14HSOUND  VELOCITY# 14 , 10. , 0. 0»CR(NU*V  +3)  , 
lCR(NUMV+4) ,10.0  ) 

CALL  AXIS  (0.0,0.0,UHWATeR  DEPTh»U#10.0,90.0»ZR(NUMV  +3)  , 

1ZR(NuMV+4) ,10.0  ) 

CALL  PLOT  (15, 0,0, 0,-3) 

ZR (NUMV  +1)  =  200 , 0*FSC 

ZR(NUMV  +2)  =  -20.0*FSC 

CR(NUMV  +1)  rVELl 

CR (NUMV  +2)  =10,0 

CALL  LINE  (CR»ZR,NUMV»1,0,0) 

CALL  LINE  ( ZBF , ZBT , 2 , 1 #  0 , 0 ) 

CALL  SYMBOL  ( 10 ,2b, ZMB, 0 , 14, 6HB0TT0M# 0 .0,6  ) 

CALL  AXIS  (0,0, 0,0, 14HS0UN0  VELOCITY » 14, 10 ., 0. 0,CR (NUMV  +1)  , 
lCR(NUMV+2) ,10.0  ) 

CALL  AXIS  (0.0,0.0,11H*ATER  DEPTH, 11 » 10 .0 , 90 . 0 »ZR (NUMV  +1)  , 

lZR(NUMV+2) ,10.0  ) 

call  plot  i  i5.o,o.o#  -3) 

8  JP  =  3 

1B8  WRITE (4 ,9) HED 

9  FORMAT  ( lhl , 9X , 20A4 ) 

WRITE (4 , 10 ) 

10  FORMAT (1H0,15X,4UZ, FT, 5X,8HC, FT/SEC) 

11  WRlTE(4,12)ZR(l) »CR(I> 

12  FORMAT (1H0,9X#F10»1,F12,1) 

IF (DABS (Z ( 1 ) -ZM) -. 01 ) 14,14,13 

13  I  =  1+1 
JP  =  JP+1 

IF ( JP-JPM) 11,8,8 

14  READ (3# 15) N,  UM  ,FQ 

15  FORMAT (110, 2F 10, 3) 

NNN  =  N 

151  READ  (3,150)  ISFQ,  IEFU,  NMOD,  INCF 
150  FORMAT  (4110) 

IF  (IRP.NE.l)  GO  TO  120 

READ  (3,121)  IRST,  IREN,  IRic#  NPS,  ZS,  ZRC  ,FMI 

121  FORMAT  (4U0,  3F10.3) 

ZS  =  ZM  -  ZS 

ZRC  =  ZM  -  ZRC 

REaD  (3,122)  (DD(J)»  J=  l.NMOD) 

122  FORMAT  (4F10.5) 

120  DO  154  Mr  l#NMOD 

DO  155  L=  ISFO#  IEFO.-INCF 
IF  ( I VOP.Nt . 1 )  GO  TO  162 

FO  =  L 

162  P2  =  6.260185306 

W2  =  ( P2*FU) * (P2*FG ) 

N  =  NNN 
FMQD  =  m 
IRtL  =  0 
IF (N ) 29 , 4 , 29 

29  K20  =  ( (W2*(CB-C(1) )*(CB+C(1) ) ) / (C ( 1 ) *C ( 1 ) *CB*CB) ) 

DK2=  . 099*k20 
K2  =  -K20 
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31  K2  =  K2  ♦  0K2 

IF (K2-K20 ) 16,32,32 

32  WRITE(4,33) 

33  FORMAT ( 1HO, 13HN0  MODE  FOUND) 

IEFQ  =  L  ♦  1NCF 

IF  (IVOP.Eu.l)  GO  TO  160 
GO  TO  111 

lb  FLB  =  FLOAT (N) 

H  =  ZM/  DBLE(FLB) 

ZP  =  0. 

J  =  1 
I  =  1 

17  IF  (ZP-  Z(I  +  D)  18,  18,19 

19  I  =  I  +  l 
GO  TO  17 

18  CZ  =  (C(I)*(ZU+1)-ZP)+C(I+1)*<ZP-Z(I)))/<Z(I+1)-Z<I)) 

CCC(J)  =  CZ 

G  =-(W2*(C(l)-CZ)*(C(l)+CZ) )/<CZ*CZ*C(l)*C(l) ) 

F( J)  =  K2  -  G 
IF(ZP-ZM)  20,21,21 

20  J=J+1 

fla  =  float (J  -l) 

ZP  =  ZM*DBLE(FLA)/DBLE(FLB) 

GO  TO  17 

21  U(l)  =  l.o 

UP(1)  ={RO/RB)*DSQRT  ( <W2*(CB-C(1) )*<CB+C<1) ) ) /( C ( 1 ) * <C ( 1 ) *CB*CB)  ) 
1-K2) 

J  =  1 

22  JP  =  3 

2b  FLA  =  FLOAT (J  -1) 

FLB  =  FLOAT(N) 

ZP  =  ZM*l>BLE(FlA)/DBL£(FLB) 

IF ( J-N-l 127,34,27 
27  IF(DABS(U(J) )-UM)28,34,34 
2d  J  =  J+l 
JP=  JP+1 
H2=  H*H 

DEN  =  1.0  ♦  (H2*F( J) )/6.0 

U(J)=  ( (1.0-(H2/3,0)*F(j-d )*U(J-1)+H*UP(J-1) ) /DEN 
UP(J)=( ( 1 . 0- (H2*F ( J) )/3.0)*UP(J-l)-(F(O-l)+F( J) - (H2*F ( J) *F < J-l ) )/ 
16.0) *H*U ( J-l ) *  0,5) /DEN 
IF  (JP-jPM)25,22,25 

34  JM  =  J 
CALL  COUNT 

19b  IF  ( MS-M) 3i , 35,35 

35  K2  =  K2  -  DK2 
DK2=  0K2  /  10.0 

IF (DK2-  , 000001*K20  *10** (-IEX )) 36, 36, 31 
3o  J  =  2 

S  =  H*U(l)*U(l>/2.0 

S2=  H*O(l)*U(l)/(2.0*CCC(l)*CCC(l) ) 

37  IFlJ-JMS)  jti,39,38 
3d  S  =  H*u ( J ) *U ( J) +S 

S2-  H*0(J)*U(J)/(CCC(0)*CCC(J) )  +S2 

J  -  J+l 
60  TO  J7 

3^  b  =(m*0( J)*U( J)/2.0)  +  b 

So=  H*o t J i *u( J) / (2. 0*CCC ( J) +CCC ( J) )  +52 


UNCLASSIFIED 
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IF  (IVOP.EO.l)  GO  TO  164 
J  =  1 

42  JP  =  3 

WRiTe(4,23)HED'N,K2,UM,fQ,M 

23  FORMAT(1h1,5X'20A4,/,2X,3hN  =,I5,2X,4hk2  =,f15.8»2X»4HUM  =>f6.1> 
12X , 6HFREQ  =  ,  F8, 1 , 2X  >  3HM  =  ,I5> 

WRITE ( 4 , 47 )  ZM,CB,RO,RB  . 

47  FORMAT ( 3X , oHZ^AX  =  ,Fl0.3,4HCB  =  ,F10.3,4HRO  =  »Fl5.8»4HRB  =  » 
1F15.8) 

WRITE(4,24) 

24  FORMAT (1H0,17X,  4HZ ,FT '  lOX , 1HU, 10X,  5HDU/0Z, 11X,  4HF ( 2 ) ) 

43  ZP  =  ZM*(  FLOAT (J-l)/  FLOAT (N)) 

ZQ  =  ZM  -  ZP 

WRITE (4,26 )ZQ,U(J) , UP( J) , FI J) 

IF(J-N-1)45,44,45 

45  IF ( DABS ( U ( J ) ) -UM >46,44,44 

46  J  =  J+l 


JP  =  JP+1 

IF ( JP- JPM  >43,42,43 
44  WRITE ( 4, 40 )  S 

40  FORMAT  (1HO,10X,  3HS  =  ,F15.5) 

164  UMAX  =0,0  ' 

NM  =  N  +1 
DO  104  K=  1*  NM 

IF  (  DABS  (U (K ) )  -  UMAX)  104,104,103 

103  UMAX  =  DABS (U (K)  ) 

104  CONTINUE 

DO  105  K=  1,  NM 

105  UU (K)  =  U (K) /UMAX 
DO  106  K=  1,  NM 

106  ZZP(K)  =  ZM*(FLOAT(K  -l)/FLOAT(N)  ) 

GAM  =  DSQRTIK20  -  K2) 

BA  =  ( Ro*C ( 1 ) ) / (2. 0*RB*CB*GAM*S  ) 

WRITE (4 , 4i )  BA,6AM 

41  FORMAT ( 1H0 , 10X, 5HB/A  =  ,F15.5,5X,7H6AMMA  =  ,F15,5  ) 

SS=  RO*S/(UMAX*UMAX)+RO*  RO*U(1)*U(1)/(UMAX*UMAX*2,0*GAM*RB) 

SSS=  R0*S2/(UMAX*UMAX>  +  RO*RO*U(l)*U(l)/(UMAX*UMAX 

1*2 . 0*6AM*RB*CB*CB) 

IF(IVOP.EQ.l)  60  TO  153 
E  =  DEXP(-6AM*H) 

ZP  =  0.0 
V  =  1.0 
NN  =  1 

ZPP(l)  =  0.0 

VV(1)  =  UU(1) 

WRITE ( 4, 23) HED, N,K2 , UM, FQ, M 
WRITE (4,47)  ZM,CB,R0,RB 
WRITE ( 4,48 ) 

48  FORMAT(1HO,17X,4HZ,FT,10X,1HU  ) 

49  WRITE ( 4 ,26) ZP , V 

26  FORMAT ( 1H0 ,F22 , 2 , F13. 3 ,F13,3 ,F16 .5 ) 

50  V  =  V*E 

NN  =  NN  +1 
VV (NN)  =  VV  (NN  -1 ) *E 
ZPP(NN)  =  ZP  -  H 

ZP  =  ZP-H 

IF  (ZPP(NN)  50.0)  110,  110,  49 

110  CALL  PLOT  (  5. 0,0.0,  -3) 
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DO  164  NINX  =  lr  NN 
164  ZPP(NINX)  =  ZM  -  ZPP(NINX) 

NINV  =  N  +  1 
DO  183  NINW  =  If  NlNV 
163  ZZP(NINW)  =  ZM  -  ZZP(NINW) 

VV (NN  +1)  =  -1,0 

VV (NN  +2)  =  0,5 

ZPP ( NN  +1)  =  200,0*FSC 
ZPP (NN  +2)  =  -20.0*FSC 
ZZP ( N  +2)  =  200 , 0*FSC 

ZZP ( N  +3)  =  -20,0*FSC 

UU(N+2)  =  -1,0 

UU(N+3)  =  0.5 

CALL  AXIS  (O.OrO.Or  9HAMPLITU0Er  9f 4 . 0 . 0 . 0 . UU(N+2 ) r UU(N+3) , 10 . 0 ) 
CALL  AXIS  (O.OrO.Or 16HWATER  DEPTH  (FT) r 16 r 10 . 0 r 90 . 0 r ZZP (N+2) , 

1  ZZP (N  +3) rlO.O) 

CALL  LINE  (UUr  ZZPr  N+l,  l.OrO) 

CALL  LINE  (VV r  ZPPr  NNr  l.OrO  ) 

CALL  LINE  ( ZBF r  ZBT r  2 r 1 r  0  r  0 ) 

CALL  SYMBOL  (10.25rZMBrO.14.6HBOTTOMrO.Or6  ) 

CALL  SYMBOL (1.3r9.50r0,14r  9HAMPL ITUDE  r  0 . 0  r  9 ) 

CALL  S YMBOL ( 1 . 7  r  9. 25 r  0 , 1 0  r  faHVERSUS  r  0 . 0  r  6 ) 

CALL  SYMBOl (1.6r9,00r0, 14r 5HDEPTHr 0 ,0 r 5) 

CAlL  SYMBOL (0.6r8,75r0 •  14  r 19HFREQUENCY  HZr0.0rl9) 

CALL  NUMBER ( 2 . 2  r  8 . 75  r  0 , 14  r  FQ  r  0 . 0  r -1 ) 

CAlL  SYMBOL ( 1 . 4 r  8 . SO  r  0 , 14  r  4HM00E  r  0 . 0  r  4 ) 

CALL  NUMBER (2.1r8,S0r0, 14  r FMODr  0 . 0  r-1 ) 

NUM  =  NUM  +1 
153  XMP(l)  =  0.0 

ZMP ( 1 )  S  0.0 
NM  =  N  +1 
CALL  PLOT  (15.OrO.Or-3) 

NSTR  =  1 
NEND  =  H  +1 
IDZ  =  0 

IF (F < 1 ) .GT.0.0)  GO  TO  179 

181  DO  170  K=  NSTRrNEND 

IF(F(K) .GT.0.0)  GO  To  171 

170  IRTB  =  K 
GO  TO  109 

171  C0S2  ( IRTB  +1)  =  FURTb  +1)*CCC(IRTB  +1)*CCC(IRTB  +l)/#2 
CV  =  CCCIIRTb  +1)/DS(jrt(1.0  -  C0S2  ( IRTB  +1)) 

ZMP  ( 2 )  =  rl*  (CV  -  CLCdRTB  +1 ))/( CCC  ( IRTB)  -CCC  ( IRTB  +D) 

XMP ( 2 )  =  (ZMP(2)  -  ZMP(i) )* (CCC (IRTB  +1)  +  CV) / (CV* (OSORT (C0S2 

1  (IRTB  +1) ) ) ) 

ZMP ( 1 )  =  FLOAT(IRTo)*H  -  ZMP (2) 

ZMP ( 2 )  =  ZMP(l)  +  ZMP (2) 

NM  =  Nrf  -IRTB  +1 
DO  172  K  =  3.NM 
TH1  =90.0 
KO  =  IRTB  -1  +R 

C0S2 ( KO-1 )  =  F ( KO  -l)*CCC(KO  -l)*CCC(Ky  -1)/W2 
COS2(KO)=F(KO)*CCC(KO)*CCC(KO)/W2 

IF(CCCIKO) .GE.CV)  ZMP(K)  =  ZMP(K  -1)  +H*(CV-  CCC(KO  -1 ) ) / ( CCC ( KO ) 
1  -LCCtKO  -l) ) 

IF(CCC(KO) .GE.CV)  GO  TO  152 

XMP(K)  =  XMP ( K  -1)  +H*(CCC(K0  -1)  +  CCC ( KO ))/( CV* (OSORT ( C0S2 ( KO-1 
1))  +  DSGiRT (C0S2  (KO)  )  ) ) 
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172  ZMP(K)  =  (KO  -1 ) *H 

TH2  =  (AC0S(SQRT(C0S2<NM  +JRTB  -1) ) ) )*57.2958 
N  =  NM  -1 
GO  TO  173 

179  C0S2(1)  =  ABS(F(1))*CCC<1)*CCC<1)/W2 

CV  =  CCC ( 1 ) /DSQRT ( 1 • 0  -  C0S2( 1 )  ) 

DO  140  K=  2'NM 

C0S2(K  -1)  =  ABS(F(K  -1))*CCC<K  -1)*CCC<K  -1)/W2 
C0S2(K)  =  ABS(F(K) )*CCC(K)*CCC(K)/W2 

TH1  = ( ACOS ( SORT ( C0S2 ( 1 ) ) ) ) *57 . 2958 
TH2  =(AC0S(SQRT(C0S2(NM) ) ) )*57.2958 
IF (  CCC(K) .GE.CV)  GO  TO  142 

XMP(K)  =  XMP(K  -1)  ♦  H*(CCC(K-1)  ♦  CCC(K) )/ (CV* ( DSQRT (C0S2 (K-l ) )♦ 
1  DSQRT (C0S2(K)))) 

140  ZMP(K)  =  (K-1)*H 

173  DO  141  K=  1 ,  N 

XMP (NM  +K)  =  2.0*XMP(NM)  -  XMP(NM  -  K) 

141  ZMP ( NM  +K )  =  ZMP (NM  -K) 

GO  TO  143 

142  ZMP(K)  =  ZMP(K  -1)  ♦  H* ( CV  -  CCC (K“l ) ) / (CCC (K)  -  CCC(K  -1)) 

152  TH2  =  90.0 

XMP (K)  =  XMP  (K  -1)  ♦  (ZMP (K)  -  ZMP (K“l ) ) * (CCC (K-l )  ♦  CV)/(CV*( 

1  DSQRT (C0S2(K  -1) ) ) ) 

N  =  K  “1 
NM  =  K 
DO  144  K=  1 ,  N 

XMP (N+l+K)  =  2. 0*XMP (N+l )  -  XMP(N+1-K) 

144  ZMP (NM  +K)  =  ZMP (NM  -K) 

143  XMP ( 2*NM)  =  0.0 

182  XMP(2*NM  +1)  =  1000.0 

NIZZ  =  2*NM  -1 
DO  186  NINZ  =  1*  NIZZ 
186  ZMP(NINZ)  =  ZM  -  ZMP(NInZ) 

ZMP (2*NM)  =  200 . 0*FSC 

ZMP(2*NM  +1)  =  -20 . 0*FSC 

WRITE  (4. 145)  TH1 r  TH2  ,XMP(NM  +N) 

IF  (IVOP.EQ.l)  60  TO  157 

145  FORMAT  (3F10.3) 

CALL  LINE  (XMP.ZMP.NM  +N, 1,0,0) 

IF  (IDZ.EQ.l)  GO  TO  650 
IDZ  =  1 

CALL  LINE  (ZBF.ZBT.2,1,0.0) 

CALL  SYMBOL  ( 10 .25 .IMS , 0 . 14, 6HB0TT0M, 0 . 0 ,6  ) 

CALL  AXIS  (0,0. 0,0 » 10HRANGE  (FT) » 10 , 12. 0 » 0 . 0 , XMP <2*Nm) 

1XMP(2*NM  +1),10.0) 

CALL  AXIS  (0.0»0.0> 16H  DEPTH  (FT ) . 16 . 10 . 0 .90 . 0 . ZMP (2*NM) . 

1  2mP(2*nM  +D.10.0  ) 

CALL  SYMBOL (5.0.9,50.0,14. 14HRAY  EQUI VALENT . 0 . 0 . 14) 

CAlL  SYMBOL (5.0.9,25.0, ^4, 19HFKEQUENCY  HZ. 0.0.19) 

CALL  NUMBLK (6 .4 .9 ,^5 »0,14»FQ.0,0»-1) 

CALL  SYMBOL  (5, 5.9,00.0, 14, 4HM0DE  .0,0.4) 

CALL  NUMBER (6.2,9.00,0, 14, PMOD. D ,0.-1) 

NM  =  NtND 

650  NSTR  =  IRTB  +  NM 

ZMP ( 1 )  =0.0 

IF(NSTR.GE.NENU)  GO  to  189 
GO  TO  lfcl 

189  CALL  PLOT  ( 15. 0 » 0 . 0 ,-3) 
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IF  (IVOP.NE.l)  60  TO  154 
157  T1(L)  =  TH1 
T2(L)  =  TH2 

FMH  =  SORT (P2*FQ  /(30.4e*CV)  ) 

PM(L)  =  SORT (P2) *  R0/(FMH*2,0*(SS)  ) 

GVEL(L)  =  SS/(CV*SSS) 

FQP(L)  =  L 

WRITE  (  4,166)  PM (L)  *  FMH,  SS  ,  GVEL(L) 

166  FORMAT (4F10 #3) 

IF  (IRP.EG.l)  PM(M)  =  PM(L) 

FMK(M)  =  P2*FQ/CV 

IZS  =  ZS*FLB/ZM  +1,0 

IZR  =ZRC*FlB/ZM  +1,0 

USR(M)  =  U ( IZS) /UMAX 

URK(M)  =  U(IZR)/UMAX 

PMM(M)  =  PM(M)*USR(M)*URR(M)/RO 

155  pvehl)  =  ev 

IF  (IRP.EQ.l)  60  TO  154 

160  IDF  =  (ISFQ  -  IEFQ) /INCF  +1 
LIFO  =  IDF  -1 

DO  161  LLL=  1»  IDF 
ILl  =  (LLL  -1)*INCF 
FQP(LLL)  =  FQP ( IEFQ+lLL) 

GVEL(LLL)  =  GVELdEFQ+ILL) 

Tl(LLL)  =  TKIEFU+ILL) 

T2(LLL)  =  T2(IEF0+ILL) 

PM  (LLL)  =  PM  ( IEF.Q  +  ILL) 

161  PVlL(LLl)  =  PVEL  ( IEFQ+ILL) 

WRITE(4,lb3)  ( PVEL ( J) »  GVEL(J),  FQP(d) ,  PM(J)»  J=  1,IDF) 
163  FORMAT  (4F10.3) 

FGP(LIFG  +2)  =  0.0 

FQP(LIFQ  +3)  =  ISFQ/10 


PVEHLIFQ  +2) 
PVEHLIHS  +3) 
GVtL (LIFQ  +2) 
GVtULlFQ  +3) 
PM (LIFO  +2) 
PM (LIFO  +3) 


3500,0 
250.0 
3500.0 
250.0 
0.0 
0.1 

CALL  PLOT  (15. 0,0. 0,-3) 

CALL  LINE  (FQP.PVEl,lIFQ  +1,1,10,14) 
call  LINE  (FQP,GVEL,LlFQ  +1,1,10,28) 

CALL  LINE  (FQP,PM,lIFQ  +1,1,10,4) 

CAlL  axis  (0.0,0.0,14HFREOUENCY  (HZ) ,14,10. 0,0.0, FQP(LIFQ  +2), 

1  FGIP(LIFQ  +3)  ,10.0) 

CALL  AXIS  ( 0 . 0 , 0 , 0 , 27HSOUND  VELOCITY  (FT  PER  SEC )  ,  «i7, 10 . 0 , 90 . 0 , 

1  PVEL  (LIFQ  +2) ,PVEL(LIFQ  +3), 10.0) 

CALL  AXIS(10.0,0.0,18HEXITATIOn  PRESSURE. -18, 10. ,90. 0»PM(LIFQ  +2), 
1  PMLIFG  +3), 10.0) 

CALL  SYMBOL (2.4,9.50,0.14, 38HS0UND  VELOCITY  AND  EXCITATION  PRESSUR 
IE, 0.0, 36) 

CALL  SYMBOL (4 . 7,9.25,0.10, 6H VERSUS ,0.0,6) 

CALL  SYMBOL (4.4,9.00,0. 14.9HFREGUENC Y ,0.0,9) 

CALL  SYMBOL (4 • 5»  8.75,0.14, 4HM0DE, 0.0,4) 

CALL  NUMBlK (5.2,8.75,0.14,FMOD,0.0,-1) 

CALL  SYMBOL (4.7,8.50,0.10, 14HPHASE  VELOCITY , 0 . 0 , 14 ) 

CAlL  SYMBOL (4. 2, 8. 50, 0.10, 14,0. 0,-1) 

CALL  SYMBOL (4.7,8. 25, 0.10, 14HGR0UP  VELOC ITY , 0 . 0 , 14 ) 

CALL  SYMBOL (4. 2, 8. 25, 0.10, 28,0. 0,-1) 
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CALL  SYMBOL (4.7  *8 . 00  * 0,10* 19HEXCITATI0N  PRESSURE* 0 .0 .lg) 

CALL  SYMBOL (4. 2*8.00 *0,10 *4,0.0 *-l) 

TKLIFQ  +2)  =  90.0 

TKLIFG  +3)  =  “10.0 

T2(LIFQ  +2)  =  90.0 

T2(LIFG  +3)  =  -10.0 

CALL  PLOT  (15. 0,0.0* -3) 

CALL  LINE  (FQP*  Tl*LlFG  +1*1*10*4) 

CALL  LINE  (FQP*  T2.LIFQ  +1,1,10*5) 

CALL  AXIS  ( 0 • 0 * 0 . 0 * 14HFREQUENCY  (HZ) *-14*10. *0.0*FQP(LIFQ  +2), 

1  FQP(LIFQ  +3) *10.0) 

call  axis  (0. o*o.o* 28H angle  of  incidence  (degrees) *28*  9.0*90.0, 

1  TKLIFQ  +2)*  TKlIFQ  +3)*10.0) 

CALL  SYMBOL (3.7*9,50*0,14 » 18H ANGLE  OF  INCIDENCE* 0 . 0 , 18) 

CALL  SYMBOL (4.7*9,25*0,10* bHVERSUS* 0 ,0*6) 

CALL  SYMBOL (4.4*9,00*0,14* 9HFREQUENC Y *0,0*9) 

CALL  SYMBOL (4.5*8.75*0.14, 4HM0DE *0.0*4) 

CALL  NUMBER (5.2*8,75*0,14*F MOD  *  0 . 0  * - 1 > 

CALL  SYMBOL (4,5* 8 ,50 *0,10* 25HB0TT0M  ANGLE  OF  INCIDENCE* 0 . 0 *25) 

CAlL  SYMB0H4. 0*8. 50*0. 10*4,0. 0,-1) 

CALL  SYMBOL (4. 5* 8. 25 *0,10* 26HSURFACE  ANGLE  OF  INCIDENCE* 0 . 0 * 26) 
CALL  SYMBOL (4. 0*8. 25* 0.10* 5*0.0 *-l) 

_Sta15+  CALL  PLOT  ( 15. 0 * 0 . 0 * -3) 

IF  (IRP.NE.l)  GO  TO  111 
123  DO  124  IP=  1  * NPS 
READ  (3,125)  NMS 

125  FORMAT  (110) 

READ  (3*126)  (MD‘S(J)*  J=1*nMS) 

126  FORMAT  (6110) 

DO  131  IR  =  IRST*  IREN,  iRic 
IREO  =  (IR  -IRSTJ/IRIC  +1 
PCT(IREO)  =  0.0 

131  PST ( IREO)  =  0.0 

DO  127  IQ=  1*NMS 
M=  MDS (IQ) 

DO  129  IR  =  IRST*  IREN,  IRic 
RR  —  I R 

IREO  =  (IR  -IRST)/IRIC  +1 
DDT (M)  =  10.0** (“Do ( M) +RR/20 . 0 ) 

PST(IREO)  =  PST( 1REO)+PMMi(M) *DDT(M)  *SIN(FMR(M)*RR-P2/8,0) 

PCT ( IRtO ) S  PCT(1RE0)+PMM(M)*DDT(M)  *COS(FMR(M)*RR-P2/6,0) 

129  CONTINUE 

127  CONTINUE 

DO  132  IR  =  IRST*  IREN*  IRIC 
RR  —  IR 

IREO  =  (IK  -IRST)/IRIC  +1 

PPT(IREO)  =  SQKT (PCT ( IREO) *PCT( IREO)  +  PST ( IREO) *PST ( IREO) ) *RO/ 
1SQRT ( RR/3 . 0 ) 

FR  =  FQ/ 1000 . 0 

A  =  (  0 . 1*FR**2/ (1.0  +FR**2) )  +  (40. 0*pR* *2/(4100.0  +Fk**2)) 

1  +  ( ,00027S*FR**2) 

PPT ( IRfc.0  )  =-20.0*ALOG10(PPT(lREO))  +A*IR/3000.0 

132  RP ( IREO  )  =  FLOAT(Ik)/6000.0 

WRITE  (4*128)  (PPT ( J) *  J=  l.IREO  ) 

12b  FORMAT  (10F10.5) 

RP ( IREO  +1)  =  0.1533 

RP( IREO  +2)  =  FMI 

PPT ( IREO  +1)  =  130.0 


88 


TR  4319 


PPTdREO  +2)  =  -20.0 

ZS  =  ZM  -  ZS 
ZRC  =  ZM  -  ZRC 

CALL  LINE  (RP»PPT » IREOr 1.10.4) 

CALL  AXIS(0.0,0.0,13HRANGE  (MILES) » l3, 20. 0, 0. 0 ,RP( IREO  +1). 
1  RPdREO  4-2), 10.0) 

CALL  AXIS(0.0»0.0»21HPROPAGATION  LOSS  ( DB) . 21 , 10 . 0 , 90 . 0 , 

1  PPTdREO  +1), PPTdREO  ♦2),io.O) 

CALL  SYMBOL (8.9, 9, 50, 0.14, 16HPR0PA6AT ION  LOSS, 0.0, 16) 

CALL  SYMBOL (9.7 ,9.25,0.10, 6H VERSUS ,0.0,6) 

CALL  SYMBOL (9. 7, 9. 00, 0,14, 5HRANGE ,0.0,5) 

CALL  SYMBOL (8, 9, 8, 75,0. 14 , 18HFREQUENCY  HZ,0.0,18) 

CALL  NUMBER (10, 3,6.75, .14,FQ,0.0,-1) 

CALL  SYMBOL (8,6, 8. 50 ,0.14 , 20HSOURCE  DEPTH  FT»0,0»20) 

CALL  NUMBER (10 , 5, 6,50* 0, 14, ZS » 0 , 0 ,— 1 ) 

CALL  SYMBOL (8.6,8.25, 0. 14, 22HRECE I VER  DEPTH  FT, 0.0,22) 

CALL  NUMBER (10. 6, 8. 25, 0.14, ZRC, 0. 0,-1) 

CALL  SYMBOL  (8.6,8.0,0.14, 5HM0DES, 0. 0,5) 

DO  187  J=  1 ,NMS 
FDS(J)  =  MDS(J) 

187  CAlL  NUMBER  (9.1  +0 . 5*J, 8. o , 0 . 14, FDS ( J) , 0 . 0 ,-l ) 

CALL  PLOT  (15. 0,0. 0,-3) 

124  CONTINUE 

CALL  PLOT  (15. 0,0. 0,-3) 

111  CALL  PLOT  (15. 0,0. 0,-3) 

END 

SUBROUTINE  COUNT 
DIMENSION  U(500)‘ 

DOUBLE  PRECISION  U 

COMMON  MS, JM,U 

COMMON  UMS 

MS=0 

J=1 

IS=1 

JMb=l 

5  IF  (U(J))  1,2*5 
$  IS1=IS 
IS=1 

7  IF  (IS-lSl)  4,2,4 
4  MS— MS-f  1 
JMS=J 
2  J=J+1 

IF  (U-JM-1)  5,0,5 
1  IS1=IS 
IS=0 
GO  TO  7 
fa  RETURN 
END 
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Appendix  B 

LIST  OF  PROGRAM  S1548 


DOUBLE  PRECISION  F,U,UP 

DOUBLE  PRECISION  K2,  *2,P2,H,ZP,  CZ,G,H2,DEN 

DOUBLE  PRECISION  S, K20 ,OK2 . GAM, BA 

DOUBLE  PRECISION  E»V 

DOUBLE  PRECISION  UMaX 

DOUBLE  PRECISION  C0S2  ,  CCC  •  CV 

DIMENSION  HED(20)»Z(l00),c(100),F(500).UP(500)  »U(500)  , 

1  ZMP ( 1000 ) »  XMP (1000) »  C0S2 ( 500 ) »  CCC(500) 

1  •  ZZP  (500 )  *  ZPP(500),  VV  ( 500 )  »  DATA(1024)  .  UU(50o) 

1  ,PM(4000)  *  ICVP(200)  ,  IADS(200)  ,CC(500)»  2Z(500) 

1  •  ZR(500)  ,  CR ( 500 )  »RP(1000) ,ZBT(5)  ,  ZBF(5) 

1 , DD ( 100 ) fPPT(lOOO) , PCT ( 1 000 ) , PST (1000) ,MDS(100) ,USR(100) «URR(100> , 
1  PMM ( 100 ) , FMR( 100 )  *FDS(100) .DDS(IOO) .  DDT(IOO) 

1  * UVEL (4000 ) ,  PVEL(4000),  FQP(4000),  TK4000),  T2(4000) 

1 »FKB ( 1 0 0 ) .ERS(lOO) »6KB(100).ER(100).PMI(100)  * GMR ( 100 ) 

COMMON  MS.jMrU 
COMMON  JMS 

CALL  PLOTS  (  DATA ( 1 ) *  1024,  6  ) 

READ(3,1)HED 

1  FORMAT  (20A4) 

READ  (3,126)  NPCS,  NVPC  ,  NADS  ,NCOR»  NMOD 
READ  (3,126)  (ICVP(J),  US  1,  NVPC) 

READ  (3,126)  (IADS(J)>  Js  1,  NADS) 

READ  (3,19b)  IRIC ,  NPS,  ZS,  ZRC,  FMI  ,FSC  ,FMJ,  FMK 
196  FORMAT  (2110,  6F10.3) 

FINC  =  30* (NaDS*NMOD  +NVPC) 

CALL  PLOT  (FINC  ‘,0.0, -3) 

DO  111  L=  l.NPCS 
NMF  =  0 

DO  189  Js  i,nVPC 
IF  (ICVP(u) .EQ.L)  GO  TO  4 
169  CONTINUE 
GO  TO  14 

4  JPM  =  25 

FINC  s30*NADS*NMOU+(NVPc  -u  +1)*30 
call  Plot  (-finc,o.o,-3j 

READ  (3,112)  NUMV ,  VEL1  ,  IVPL  , IEX 

112  FORMAT  (  110, F10, 0,3110) 

NUM  s  0 
READ (3,2)  ZMM*CB 

2  FOkMAT  (4F10.3) 

ZMB  =  (200 , 0#FSC  -ZMM) / (20 . 0*FSC ) 

ZBT(l)  =  ZmM 
ZBT ( 2 )  =  ZMM 
ZBT ( 3 )  =  200 , 0*FSC 
ZBT (4)  =  -20 , 0*FSC 
ZBF ( 1 )  s  0.0 
ZBF (2)  s  10,0 
ZBP ( 3 )  s  o.O 
ZBF (4)  s  i.o 
DO  3  1=1,100 
ZZ( I)  =  0.0 

3  CC  ( I )  =  0.0 

1  =  1 

5  READ  (3.30)  ZZ(1),  CC ( I ) 

30  FORMAT ( 2F10 . 3) 

IF ( DABS ( ZZ ( I )  -  ZMM)  -  .01)  6,6,7 
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71=1+1 
GO  TO  5 
6  1  =  1 

IF  (IVPL.EQ.l)  GO  TO  8 
DO  130  K  =  lr  NUMV 
CRlNUMV  -K  +1)  =  CC (K) 

130  ZR (NUMV  -K  +1)  =  ZZ<K) 

ZR (NUMV  +1)  =  ZMM  -  O.Ol 
ZRtNUMV  +2)  =  ZMM  +  50,0 
CR ( NUMV  +1)  =  CB 

CR(NUMV  +2)  =  CB 

DO  185  NINY  =  1,  NUMV 
185  ZR (NINY )  =  ZMM  -  ZR(NlNY) 

ZR(NUMV  +3)  =  200.0*FSC 

ZR(NUMV  +4)  =  -20.0*FSC 

CR(NUMV  +3)  =  VEL1 

CR (NUMV  +4)  =  50,0 

CALL  LINE  (CRrZRrNUMV  +2,  1,0,0  ) 

CALL  LINE  (ZBF,ZBT,2,1,0,0) 

CALL  SYMBOL  ( 10,25r ZMB, 0 , 14 » 6HB0TT0M, 0 .0*6  ) 

CALL  AXIS  (0.0, 0.0, 14HS0UND  VELOCITY, 14, 10 . , 0 . 0 ,CR ( NUMv  +3)  , 
lCR(NUMV+4) ,10.0  ) 

CALL  AXIS  (0.0,0.0,llHWATER  DEPTH, 11  *  10 . 0 . 90 . 0 .ZR (NUMV  +3)  , 

lZR(NUMV+4) ,10.0  ) 

CALL  PLOT  (15. 0,0. 0,-3) 

ZR (NUMV  +1)  =  200,0*FSC 

ZR (NUMV  +2)  =  -20.0*FSC 

CR (NUMV  +1)  =VEL1 

CR (NUMV  +2)  =10.0 

CALL  LINE  (CRrZR, NUMV, 1,0,0) 

CALL  LINE  ( ZBF , ZBT , 2 , 1 , 0 , 0 ) 

CALL  SYMBOL  ( 10 .25 ,ZMB, 0 . 14 , 6HB0TT0M, 0 . 0 , 6  ) 

CALL  AXIS  (0.0, 0.0, 14HS0UND  VELOCITY, 14 , 10 ., 0. 0 ,CR (NUMV  +1)  , 
lCR(NUMV+2) ,10.0  ) 

CALL  AXIS  (0.0,0.0,11HWATER  DEPTH, 11 , 10 .0 , 90 .0 , ZR (NUMV  +1)  , 

1ZRINUMV+2) .10.0  ) 

CALL  PLOT  (-15. ,0.0, -3) 

CALL  PLOT  (FINC  ,0.0, -3) 

8  JP  =  3 

188  WRITE ( 4 ,9) HED 

9  FORMAT  (1H1»9X,20A4) 

WRITE (4 , 10 ) 

10  FORMAT ( 1H0 , 15X , 4HZ , FT , 5X , 8HC , FT/SEC ) 

11  WRITE(4,12)ZR(I),CK(I> 

12  FORMAT (1H0,9X*F10.1*F12,1) 

IF (DABS(ZZ( I )  -ZMM)  -  0.01)  14,14,13 

13  I  =  1+1 
JP  =  JP+1 
IF(JP-JPM)11,8,8 

14  IF  (L.NE.l)  GO  TO  197 
READ  (3,15)  N,  UM,  F® 

15  FORMAT (110, 2F 10, 3.110) 

NNN  =  N 

IF  (L.NE.l)  GO  TO  197 
DO  204  J=  1,  NM 
FKb(M)  =  0.0 
DDb(M)  =0.0 
204  ERS(M)  =  0.0 
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197 

195 


122 


200 


201 

202 


205 


205 

120 

152 


READ  (3* 195)  ZM*  1RST  ,  IREN 
FORMAT  (F10 . 0  *  2110) 

READ  (3*2)  CB  r  RO  r  RB  *XI 
READ  (3,122)  (DO(J)*  J=  1,NM0D) 

FORMAT  (4F10.5) 

IPOS  =  0 

IF  (ZRC.GT.ZM)  IPOS  =  1 
ZS  =  ZM  -  ZS 
ZRC  =  ZM  -  ZRC 

ZMB  =  ( 200 . 0*FSC  -ZM  )/(20.0*FSC) 

DO  200  J=  1 * NUMV 
C(J)  =  CC(J) 

Z(J)  =  ZZ(J)  -  ( ZMM  -ZM) 

DO  201  J=  1*NUMV 

IF (Z( J) .LT.0.1)  GO  TO  201 

GO  TO  2 02 

CONTINUE 

IF  (J.EG.l)  GO  TO  120 

C(J  -1)  =  C(J  )  +(C(J)  -C(J  -1) )*-Z(J)/(Z(J) 

Z(J  -1)  =  0.0 

JUO  =  NUMV  -J  +2 

DO  203  JJ=  1*JJJ 

C(UJ)  =  C(JJ  +J  -2) 

Z(bJ)  =  Z(JJ  +J  -2) 

JJu=  JJJ  +1 

DO  205  JK.  =  JJJ * NUMV 

C (oK)  -  0.0 

Z(JK)  =  0.0 

DO  154  Ms  1 * NMOD 

P2  =  6.283105306 

W2  =  (P2*FG)*(P2*Fu) 

N  =  NNN 


Z(J  -D) 


FMuD  =  M 

IF (N ) 29 , 4 » 29 

29  K20  =  ( (W2* (CB-C (1) ) *(CB+C ( 1 ) ) ) / (C ( 1 ) *C ( 1 ) *CB*CB) ) 
DK2=  . 099*K20 
K2  =  -K20 

31  K2  =  K2  +  uK2 

IF (K2-K20 ) 16*  32  *  32 

32  WRITE ( H  r  33) 

35  FORMAT (1H0.13HNO  MODE  FOUND) 

NMF  =  1 

DO  295  J=  1,  NMOD 

29b  DDS(J)  s  uuS(J)  +Du ( J) *FLOAT ( IRfcN) 

GO  TO  123 

Id  FLb  =  FLOAT  (N) 

H  =  ZM/  DbLE(FLB) 

ZP  =  0. 

J  =  1 
I  =  1 


17 

IF 

(ZP-  Z ( 1+1 ) )  16 

,18*19 

19 

I 

=  I  +1 

GO 

TO  17 

ia 

cz 

=  (C(I)*(Z(I+1) 

-ZP)+C(I+1)*(ZP-Z(I) 

) )/(Z(I 

ccc 

(J)  =  CZ 

G  = 

- ( W2* ( C ( 1 ) -CZ) * 

(C(1)+CZ) )/(CZ*CZ*C( 

1)*C(1) 

F ( J)  =  K2  -  G 
IF(ZP-ZM)  20,21,21 
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20  d=0+l 

FLA  =  FLOAT (d  -1) 

ZP  =  ZM*DBLE (FLA ) /DBLE ( FLB) 

GO  TO  17 

21  U(l)  =  1.0 

UP(1)  =(RO/RB>*DSQRT  ( (W2*(CB-C(1) )*(CB+C(1) ) )/(C(l)*(C(l)*CB*CB) ) 
1-K2) 
d  =  1 

22  dP  =  3 

25  FLA  =  FLOAT ( J  -1) 

FLB  =  FLOAT (N) 

ZP  =  ZM*DBLE(FLA)/DBLE(FLB) 

IF(d-N-l)27»34,27 

27  IF(DABS(U(J) )-UM)26,34,34 

28  J  =  J+l 
JP=  dP+1 
H2=  H*H 

DEN  =  1.0  +  (H2*F(d) >/6,0 

U(d)  =  ( (1.0-(H2/3.0)*F(d-l) )*U(j-l)+H*UP(d-l) )/DEN 
UP(d)=( (1.0-(H2*F(J) )/3.0)*UP(J-l)-(F(J-l)+F(d)-(H2*F( J)*F(d-l) )/ 
16.0) *H*U (d-1 ) *  0.5)/DEN 
IF  (dP-dPM)25»22,2b 

34  JM  =  J 
CALL  COUNT 

IF  (MS-M)31»35,35 

35  DO  289  J=  1 .N  ‘ 

IF  (F(J).GT.O.O)  GO  TO  290 

289  CONTINUE 

290  IJ1  =  d  +1 

DO  291  d=  IJl*  N 
IF(Ftd).LT.O.O)  GO  TO  292 

291  CONTINUE 
GO  TO  294 

292  IJ2  =  J  +1 

DO  293  J=  I j2»  N 
IF (F(J) ,GT .0,0)  GO  TO  31 

293  CONTINUE 

294  K2  =  K2  -  DK2 
DK2=  DK2  /  10.0 

IF(DK2-  ,000001*K20*10**(-IEX) )  36»36»31 

36  J  =  2 

IF  (XI.GT.0.0)  GO  TO  502 
GO  TO  501 

502  READ  (3,500)  (U(Ii)  •  11=  1,N) 

500  FORMAT  (8D10.5) 

501  S  =  H*U(l)*U(l)/2.0 

S2=  H*J(1)*U(1)/(2.0*CCC(1)*CCC(1) ) 

37  IF (J-JMS >36,39, 38 

38  S  =  H*d(d)*U(d)+S 

S2=  H*U(J)*U(J)/(CCC  (>J)  *CCC  ( J)  )  +S2 

d  =  d+1 
GO  TO  57 

39  S  =(H*U(d)*U(d)/2.0)  +  S 

S2=  H*U(d)*U(d)/(2.0*CCC(d)*CCC(d) )  +S2 

d  =  1 
42  dP  =  3 

DO  206  KK=  1,  NADS 
IF(IAOS(KK) .EQ.L)  GO  TO  208 
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206  CONTINUE 
GO  TO  43 

206  WRITE(4,23)HED,N,K2,UM,FQ,M 

23  FORMAT (lHl»5X»20AW»2Xr3HN  =rl5r2Xr4HK2  =,Fl5.8,2X,4HUM  =*F6.1. 
12X.6HFREQ  =  »F8.1»2X»3HM  =rl5) 

WRITE(4,47)  ZMfCBfROfRB 

47  FORMAT (3X»6HZMAX  =  »F10.3r4HCB  =  »F10.3»4HRO  =  *Fl5.8»4HRB  =  * 
1F15.8) 

WRITE(4,24) 

24  FORMAT ( 1H0 , 17X ► 4HZ »FT » lOXr 1HU» 10X » 5HDU/DZ» 11X» 4HF( Z) ) 

43  ZP  =  ZM*(  FLOAT(J-l)/  FLOAT(N)) 

ZQ  =  ZM  -  ZP 

DO  207  KK=  1,  NADS 

IF ( IADS (KK) . EQ.L)  GO  TO  209 

207  CONTINUE 
GO  TO  44 

209  WRITE  (4. 26)  ZQ» U  ( J) » UP  ( J)  ,F(  J) 

IF ( J-N-l ) 45»  44 » 45 

45  IF (DABS (U ( J) ) ”UM) 46»  44 >44 

46  J  =  J+l 
JP  =  JP+1 

IF ( JP-JPM) 43*42 » 43 

44  WRITE (4  .  40  )  S 

40  FORMAT  (lHO.lOXr  3HS  =  ,F15.5) 

164  UMAX  =0.0 

NM  =  N  +1 
DO  104  K=  1.  NM 

IF  (  DABS  (U(K))  -  UMAX)  104,104,103 

103  UMAX  =  DABS ( U ( K )  ) 

104  CONTINUE 

DO  105  K=  lr  NM 

105  UU (K )  =  U (K) /UMAX 
DO  106  K=  1,  NM 

106  ZZP(K)  =  ZM* (FLOAT (K  -1)/FL0AT(N)  ) 

GAM  =  USQRT (K20  -  K2) 

BA  =  (RO*C (1) )/(2.0*RB*CB*GAM*S  ) 

'~<\ITE(4»41 )  BA r GAM 

41  FORMAT ( 1H0» 10X» 5HB/A  =  »F15. 5» 5Xr 7HGAMMA  =  ,F15.5  ) 

S»S=  R0*S/(UMAX*UMAX)4R0*  RO*U(1)*U(1)/(UMAX*UMAX*2.0*GAM*RB) 

OSS=  R0*S2/(UMAX*UMAX)  ♦  RO*RO*U (1)*U(1) / (UMAX* UMAX 

l*.2.0*GAM*Rb*CB*CB) 

E  =  DEXp(-faAM*H) 

Z P  =0.0 

V  =  1.0 
NN  =  1 

ZPP(l)  =  0.0 

VV(1)  =  UU ( 1 ) 

WRITE(4.23)HED,NrK2rUM»FQ,M 
WRITE ( 4 » 47)  ZM»Cb*ROrRB 
212  DO  210  Kk=  1»NADS 

IF( IAQS(KR) .EQ.L)  GO  TO  211 

210  CONTINUE 
GO  TO  bo 

211  WRITE ( 4 . 48 ) 

46  FORMAT ( 1H0  » 17X  »  4HZ »  FT  * 10X  » 1HU  ) 

49  WRITE(4<26) ZP*  V 

26  FOKMAT(1HO,F22.2»F13.3,F13.3,F16.5) 

50  V  =  V*E 
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NN  =  NN  +1 
VV(NN)  =  WINN  -1)*E 
ZPP(NN)  =  ZP  -  H 

ZP  =  ZP-H 

IF  (ZPP(NN)  +  50.0)  110'  110'  212 

110  DO  190  J=  1'  NADS 

IF  (  IADS(J).EG.L)  GO  TO  193 
190  CONTINUE 
GO  TO  192 

193  FINC  =  30* (NaDS*NMOD  -J  +1)  -(NMOD  -M)*30 
CALL  PLOT  (-FlNC'0.0»-3) 

DO  184  NINX  =  1'  NN 
184  ZPP(NINX)  =  ZM  -  ZPP(NINX) 

NINV  =  N  +  1 
DO  183  NINW  =  1'  NINV 
183  ZZP(NINW)  =  ZM  -  ZZP(NINW) 

WINN  +1)  =  -1,0 

WINN  +2)  =  0.5 

ZPPINN  +1)  -S  200.0*FSC 
ZPPtNN  +2)  =  -20,0*FSC 
ZZPIN  +2)  =  200.0*FSC 

ZZPIN  +3)  =  -20 . 0*FSC 

ZBT(l)  =  ZM 
ZBT12)  =  ZM 
UUIN+2)  =  -1.0 

UU ( N+3)  =  0.5 

CALL  AXIS  tO.0'0.0.  9HAMPLITUDE'  9» 4 . 0 » 0 . 0 » UU (N+2) ,UU (N+3) , 10 .0 ) 
CALL  AXIS  (0.0'0.0'16HWATER  DEPTH  (FT) » 16. 10 .0 ' 90 .0 ,ZZp (N+2) , 

1  ZZP ( N  +3) '10.0) 

CALL  LINE  (UU.  ZZP'  N+l '  l.O'O) 

CALL  LINE  (VV»  ZPP.  NN,  1,0.0  ) 

CALL  LINE  (ZBF'ZBT'2'1'0'0) 

call  SYMBOL  ( 10 ,25. ZMB , o . 14 .6HB0TT0M. 0 ,0.6  ) 

CALL  SYMBOL ( 1 . 3 . 9 . 50,0.14, 9H AMPLI TUDE ,0.0,9) 

CALL  SYMBOL (1.7,9.25,0 • 10 '6HVERSUS, 0 .0 , 6) 

CALL  SYMBOL ( 1.6,9.00,0 • 14, 5HDEPTH, 0 .0,5) 

CALL  SYMBOL (0 .8 , 8, 75, 0 , 14 , 19HFREQUENCY  HZ, 0.0.19) 

CALL  NUMBER (2.2, 8, 75, 0 .14,FQ.0.0,-1) 

CALL  SYMBOL (1.4,8.50,0 • 14 ,4HM0DE ,0.0,4) 

CAlL  NUMBER (2. 1,8. 50,0. 14, FMOD,0. 0,-1) 

192  NUM  =  NUM  +1 
153  XMP(l)  =  0.0 

ZMP ( 1 )  =  0.0 

NM  =  N  +1 
CALL  PLOT  (15. 0,0. 0.-3) 

NSTR  =  1 
NEND  =  N  +1 
IDZ  =  0 

IF (F(l) .GT.0.0)  GO  TO  179 
681  DO  170  K=  NSTR,NEND 

IF(F(K) .G7.0.0)  GO  TO  171 

170  IRTB  =  K 
GO  TO  bb9 

171  C0S2  ( IHTB  +1)  =  F  ( IRTB  +1)*CCC(IRTB  +1)*CCC(IRTB  +l)/*2 
CV  =  CCC ( IRTB  +1)/DSgrt(1.0  -  C0S2(IRTB  +1)) 

ZMP (2 )  =.  H*(CV  -  CCC (IRTB  +1 ))/ (CCC ( IRTB) -CCC ( IRTB  +1)) 

XMP (2)  =  (ZMP(2)  -  ZMP (1) )*(CCC(IRTB  +1)  +  CV) / (C V* ( DSGRT (C0S2 

1  (IRTB  +1)))) 
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ZMP<1)  =  FLOAT(IRTB)*H  -  ZMP ( 2 ) 

ZMP (2 )  =  ZMP(l)  +  ZMP ( 2 ) 

NM  =  NM  -IRTb  +1 
DO  172  K  =  3.  NM 
THl  =90.0 
KO  =  IRTB  -1  +K 

C0S2 ( KO-i )  =  F ( KO  -1 ) *CCC ( KO  -l)*CCC(KO  -1J/W2 
C0b2(K0)=F(K0)*CCC(K0)*CCC(K0)/W2 

IF (CCC (KO) .GE.CV)  ZMP(K)  =  ZMP(K  -1)  +H*(CV-  CCC(KO  -1 ) ) / (CCC (Ko) 
1  -CCC ( KO  -1)  ) 

IF(CCC(KO) .GE.CV)  SO  TO  152 

XMP(K)  =  XMP(K  -1)  +H* (CCC (KO  -1)  +  CCC (KO ) ) / (CV* (DSORT (C0S2 (KO-1 
1) )  +  DSGRT (C0S2(K0) ) ) ) 

172  ZMP(K)  =  (KO  -1 ) *H 

TH2  =  ( ACOS ( SORT (CoS2 (NM  +IRTB  -1 ) ) ) ) *57. 2958 
N  =  NM  -1 
GO  TO  173 

179  C0S2 ( 1 )  =  ABS (F ( 1 ) ) *CCC ( 1 ) *CCC ( 1 ) /W2 

CV  =  CCC ( 1 ) /DSQRT ( 1 • 0  -  C0S2(1)  ) 

DO  140  K=  2»NM 

COS2 (K  -1)  =  ABS(F (K  “1))*CCC(K  -1)*CCC(K  -1)/W2 
COS2  ( K )  =  ABS(F(K)  )*CCC(M*CCC(K)/W2 

THl  =(ACOb(SQRT <C0b2(l) ) ) )*57,2958 
TH2  =(AC0b(SQRT(C0S2(NM) ) ) >*57.2958 
IF (  CCC(K) .GE.CV)  GO  To  142 

XMP(K)  =  XMP(K  -1)  ♦  H* (CCC ( K”1 )  ♦  CCC ( K) ) / (CV* (0SORT ( C0S2 (K-l ) ) + 

1  DSQRT(C0S2(K) ) ) ) 

140  ZMP (K )  =  (K-1)*H 

175  DO  141  K=  1»N 

XMPtNM  +K)  r  2.0*Xf.P(NM)  -  XMP(NM  -  K) 

141  ZMP ( NM  +K)  =  ZMP ( NM  -K) 

GO  TO  143 

142  ZMP(K)  =  ZMP(K  -1)  H*  (CV  -  CCC  ( K-l )  )  /  (CCC  (K )  -  CCC(K  -1)) 

152  TH2  =  90.0 

XMP(K)  =  XMP(K  -1)  +  (ZMP(K)  -  ZMP (K”1))*(CCC(K-1)  +  CV)/(CV*( 

1  DSORT (COS2 ( K  -1)))) 

N  =  K  -1 
NM  =  K 
DO  144  K=  l.N 

XMP  ( N-fl-fK )  r  2,0*XNiP(N+l)  -  xMP(N-fl-K) 

144  ZMP (N  +1  +K)  =  ZMP (N  +  1  -K) 

143  XMP  ( 2*Nf.  )  -  0.0 

162  XMP (2*NP  +1)  =.  1000.0 

NI2Z  =  2*NM  -1 
DO  186  TilNZ  =  1,  NiZZ 
16b  ZMP(NINZ)  =  ZM  -  Zi-iP(MNZ) 

ZMP  ( 2*I'<M)  =  200 .  U*FSC 

ZMP ( 2*NM  +1)  =  ”20 . 0*F  SC 

WRITE  (4.145)  THl. 1H2»XMP(NM  +N) 

14b  FOkMAT  ( 3FiO . 3 ) 

DO  691  J=  i,  NADS 
IF  (  IAuS(J).EQ.U)  GO  TO  194 
691  CONTINUE 
GO  TO  157 

194  CAlL  LINE  (XMP.ZMP.NM  +N. 1,0.0) 

IF  (ID2.EG.1)  GO  TO  650 
IDZ  =  1 

CAi_L  LINE  ( Z9F  , ZBT » 2 » 1 » 0  »  0 ) 
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CAi_L  SYMBOL  (10,25,ZMB, 0 , 14 , 6HB0TT0Mf 0 . 0 » 6  ) 

CALL  AXIS  (O.Of O.Of 10HRANGE  (FT) , 10 » 12. 0 , 0 . 0 » XMP (2*NM)  * 

1XMP(2*NM  +1) ,10.0) 

CALL  AXIS  (0.0, 0.0, 16H  DEPTH  (FT ) , 16, 10 . 0 ,90 . 0 , ZMp (2*NM) , 

1  ZMP(2+NM  +1) ,10.0  ) 

CALL  SYMBOL(5.0,9.&0,0.14,14HRAY  EQUI VALENT, 0 . 0 , 14 ) 

CALL  SYMBOL (5.0,9.25,0.14, 19HFKEQUENCY  HZfO.0,19) 

CALL  NUMBER (6. 4, 9. 25,0. 14, FQ,0. 0,-1) 

CALL  SYMBOL (5.5,9.00,0.14, 4HM0DE ,0.0,4) 

CALL  NUMBER (6.2»9.00,0 • 14 ,FMOD, 0 . 0 , -1 ) 

NM  =  NEND 

650  NSTR  =  IRTB  +  NM  -i 

ZMp(l)  =  o.O 

IF(NSTR.GE.NEND)  GO  TO  689 
GO  TO  661 

689  CALL  PLOT  ( -l5. 0 , 0 . 0 ,-3 ) 

CALL  PLOT  (FINC  , 0.0> “3) 

157  T1(L)  =  TH1 
T2(L)  =  TH2 

FMH  =  SGRT (P2*FQ  /(30.48*CV)  ) 

PM (L)  =  SORT (P2)*  R0/(FMH*2.0*(SS)  ) 

GVEL(L)  =  SS/(CV*SSS) 

FQPIL)  =  L 

WRITE  (  4, 166 )  PM(L)  f  FMH,  SS  ,  GVEL(L) 

166  FORMAT (4F10.3) 

PM(M)  =  PM(L) 

FMR(M)  =  P2*FQ/CV 

FKb(M)  =  FKB(M)  +  RMR(M)*FLOAT(IREN) 

DDS(M)  =  OUS(M)  +DU(M)*FLOAT ( IREN) 

IF  (L.EG.NCOR)PMI (M)=  PM(L) 

IZR  =ZKC*FlB/ZM  +1.0 
IZS  =  ZS*FLB/ZM  +1.0 
IF(L.EG.NCOR)USR(M)=  U(IZS)/UMAX 
IF ( IPOS.EU. 1 )  GO  TO  165 
URR(M)  =  o(IZR)/UMmX 

PMM(M)  =  S^RT(PM(M)*PMI (M) )*USR(M)«URR(M)/RO 
WRITE  ( 4 , lo6 )  PMM(M) fPMI (M) ,  USR(M)f  URR(M) 

WRITE (4, 12b)  L,  NCOR 
155  PVc.HL)  =  CV 
154  CONTINUE 
123  ZS  =  ZM  -  ZS 
ZRC  =  Zm  -  ZRC 
IRtN  =  IRS1  +  IREN 
DO  124  IP=  1 , NPS 
IF  (L.NE.l)  GO  TO  198 
READ  (3,125)  NMS 

125  FOrcMAT  (110) 

READ  (3,12b)  (MDS(J)f  J=1,NMS) 

126  FORMAT  (6110) 

198  DO  131  IR  =  IRST,  IREN,  IRIC 
IREO  =  (IR  -IRST) /IRlC  +1 
PCT(IREO)  =  0.0 

131  PS1 ( IREO)  =  0.0 

IF  (NMF.Efcl.l)  GO  TO  111 
DO  127  Io=  1 f NMS 
M—  MDS(IU) 

DO  129  IR  =  IRST,  IRENf  IRIC 
IF  (L.Eo.NMOD)RI=  FLOAT(IRST) 
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GKB(M)  =  FKB(M)  -  FLOATUREN  -  IR)*FMR(M) 

DDT (M)  =  DDS(M)  -FLOAT(IREN  -IR)*DD(M) 

IF  (GKB(M)  .LE.O.O)  GO  TO  214 
GKB(M)  =  GKB(M)/(FLOAT(IR)  -Rl) 

214  RR  =  IR 

IF(L.EO.NMOD)GKB(M)=  FMR(M) 

GMR(M)  =  GKB(M)  +  ERS(M)/RR 
ER(M)  =  FMR(M)  -  GRB(M) 

IF ( IREN  -IR.GE.IRIC)  WGT  =  IRIC 

IF (IREN  -IR.LT.IRIC)  WGT  =  IREN  -IR 

ERS(M)  =  ERS(M)  +  ER(M)*WGT 

IREO  =  (IR  -IRSTJ/IRIC  +1 

IF ( IPOS, EG. 1 )  GO  TO  129 

DDKH)  =  10.0**(-DDT(M)/20.0) 

PST ( IREO)  =  PST(IREO)+PMM(M)*  DDT(M)  *SIN (GMR ( M) *RR-P2/6. 0 ) 

PCT ( IREO) =  PCT(IREO)+PHM(M)*  DDT(M)  *COS (GMR (M) *RR-P2/6. 0 ) 

129  CONTINUE 

127  CONTINUE 

IF ( IPOS . Eu. 1 )  GO  TO  113 
DO  132  IR  =  iRSTr  IREN,  iric 
RR  —  l  R 

IREO  =  (IR  -IRST ) /IRIC  +1 

PPT(IREO)  z  SQRT(PCT(lREO)*PCT(lREO)  +  PST( IREO) *PST ( IREO) ) *RO/ 
1SQRT  (RR/3 . 0 ) 

FR  =  FU/1000 . 0 

A  =  (  0.1*FR**2/(1.0  +FR**2 ) )  +  (40.0*FR **2/(4100.0  *FR**2) ) 

1  +  ( ,00027S*FR**2) 

PPTdRtO  )  =-20.0*ALOGlO(PPT(  IREO)  )  +A*IR/3000.0 

132  RP ( IREO  )  =  FLOAT(IR)/6000.0 

WRITE  (4.128)  (PPT ( J) *  J=  1 . IREO  ) 

128  FORMAT  (10F10.5) 

113  RP ( IREO  +1)  =  0.0 

RP ( IREO  +2)  =  FMI 

PPT (IREO  +1)  =  FMJ 

PPTUREO  +2)  =  FMK 

IF ( IPOS.EU. 1 )  GO  TO  114 
CAlL  LINE  (RP.PPT. IREO. 1,10,4) 

114  IF  (L.NE.NCOR)  GO  10  111 

CAlL  AXiS(0.0»0.0,13HRanGE  (MILES) ,13, 20. 0,0.0, RPdREO  +1), 

1  RPdREO  +2)  ,10.0) 

CAlL  AXIS(0.0*0.0,21HPROPAGATION  LOSS  (DB) , 21 , 10 . 0 f 90 . 0 , 

1  PPT (IREO  +1 ) *  PPT ( IREO  +2). 10.0) 

CALL  SYn.BOl (8 . 9, 9, 50 , 0 . m , 16HPR0PAGATI0N  LOSS. 0.0,16) 

CAlL  SYFBOl (9. 7, 9, 25, 0.10, 6H VERSUS, 0.0,6) 

CALL  SYMBOL (9.7,9. 00,0.14, 5HRANGE ,0.0,5) 

CALL  SYMBOL (8. 9, 8, 75, 0,14, 18HFREQUENCY  HZ,0.0,i8) 

CAlL  NUMBER (10.3,8.75*. 14,FQ,0.0,-1 ) 

CAlL  SYMBOL (8.6,8.50,0.14, 20HSOURCE  DEPTH  FT, 0.0,20) 

CAlL  NUMBLK(l0.5,6.50r0.14,ZS*0.0,-l) 

CAlL  SYNBOl (6.8,8.25,0.14, 22HRECEI VER  DEPTH  FT, 0.0, 22) 

CAlL  NUmBER(10.6,8.25*0.14,ZRC,0.0,-1) 

CAlL  SYMBOL  (8.6,8.0,0.14, 5HMODES, 0.0,5) 

DO  187  J=  1 , NMS 
FDS(J)  =  MDS(J) 

187  CAlL  NUMBER  (0.1  +0. S*j, 8.0,0. 14, FDS(J) , 0 . U , -1 ) 

199  FINC  =  1P*25 

CAlL  PLOT  (20.0,0.0,-J) 

124  CONTINUE 
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DO  220  J=  l.NPS 
220  CALL  PLOT  (-20 . 0 . 0 . 0 .-3) 
111  CONTINUE 

CALL  PLOT  (15.0  .0.0.-3) 
END 

SUBROUTINE  COUNT 
DIMENSION  U ( 500 ) 

DOUBLE  PRECISION  U 

COMMON  MS.JM.U 

COMMON  JMS 

MS=0 

J=1 

IS=1 

JMS=1 

5  IF  (U(J))  1.2.3 

3  IS1=IS 
IS=1 

7  IF  (IS-ISl)  4.2.4 

4  MS=MS+1 
JMS=J 

2  J=J+1 

IF  (J-JM-i)  5.6.5 
1  IS1=IS 
IS=0 
GO  TO  7 

6  RETURN 
END 
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